如何找到矩阵与循环收敛的时间 [英] How to find when a matrix converges with a loop
问题描述
我得到了一个矩阵:
P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE)
使用功能mpow
,rows_equal
,matrices_equal
.我想找到P^n
何时收敛,换句话说n是什么,矩阵中的所有行均相等时以及P^n = P^(n+1)
何时.
Using the functions, mpow
, rows_equal
, matrices_equal
. I want to find when P^n
converges, in other words what n is, when all the rows are equal in the matrix and when P^n = P^(n+1)
.
仅查看函数i
即可推断出n=19-21
周围的矩阵将收敛.
By just looking at the functions i
have managed to deduce that around n=19-21
the matrix will converge.
尽管,我想使用循环找到合适的n.下面是功能mpow
,rows_equal
和matrices_equal
.我知道它们的书写方式可能会有所不同,但请保持原样.
Although, I want to find the right n using a loop. Here under are the functions mpow
, rows_equal
and matrices_equal
. I know they can be written differently but please keep them as they are.
mpow <- function(P, n, d=4) {
if (n == 0) diag(nrow(P)))
else if (n== 1) P
else P %*% mpow(P, n - 1))
}
rows_equal <- function(P, d = 4) {
P_new <- trunc(P * 10^d)
for (k in 2:nrow(P_new)) {
if (!all(P_new[1, ] == P_new[k, ])) {
return(FALSE)}
}
return(TRUE)
}
matrices_equal <- function(A, B, d = 4) {
A_new <- trunc(A * 10^d)
B_new <-trunc(B * 10^d)
if (all(A_new == B_new)) TRUE else FALSE
}
现在,要编写循环,我们应该按照以下方式进行:
Now, to write the loop, we should do it something along the lines of:
首先创建一个像这样的函数:
First creating a function like so:
when_converged <- function(P) {...}
和
for (n in 1:50)
在t.ex n = 50时尝试.
To try for when t.ex n = 50.
尽管我不知道如何正确编写代码,但是有人可以帮助我吗?
Although i don't know how to write the code correctly to do so, can anyone help me with that?
感谢您阅读我的问题.
推荐答案
实际上,更好的方法是这样做:
Actually, a much better way is to do this:
## transition probability matrix
P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE)
## a function to find stationary distribution
stydis <- function(P, tol = 1e-16) {
n <- 1; e <- 1
P0 <- P ## transition matrix P0
while(e > tol) {
P <- P %*% P0 ## resulting matrix P
e <- max(abs(sweep(P, 2, colMeans(P))))
n <- n + 1
}
cat(paste("convergence after",n,"steps\n"))
P[1, ]
}
然后,当您调用函数时:
Then when you call the function:
stydis(P)
# convergence after 71 steps
# [1] 0.002590674 0.025906736 0.116580311 0.310880829 0.272020725 0.272020725
函数stydis
本质上是连续执行的:
The function stydis
, essentially continuously does:
P <- P %*% P0
直到P
的收敛.收敛性由差异矩阵的L1范数确定:
until convergence of P
is reached. Convergence is numerically determined by the L1 norm of discrepancy matrix:
sweep(P, 2, colMeans(P))
L1范数是所有矩阵元素的最大绝对值.当L1范数降到1e-16
以下时,就会发生收敛.
The L1 norm is the maximum, absolute value of all matrix elements. When the L1 norm drops below 1e-16
, convergence occurs.
如您所见,收敛需要71个步骤.现在,我们可以通过控制tol
(公差)来获得更快的收敛":
As you can see, convergence takes 71 steps. Now, we can obtain faster "convergence" by controlling tol
(tolerance):
stydis(P, tol = 1e-4)
# convergence after 17 steps
# [1] 0.002589361 0.025898057 0.116564506 0.310881819 0.272068444 0.271997814
但是,如果您检查:
mpow(P, 17)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 0.002589361 0.02589806 0.1165645 0.3108818 0.2720684 0.2719978
# [2,] 0.002589415 0.02589722 0.1165599 0.3108747 0.2720749 0.2720039
# [3,] 0.002589738 0.02589714 0.1165539 0.3108615 0.2720788 0.2720189
# [4,] 0.002590797 0.02590083 0.1165520 0.3108412 0.2720638 0.2720515
# [5,] 0.002592925 0.02592074 0.1166035 0.3108739 0.2719451 0.2720638
# [6,] 0.002588814 0.02590459 0.1166029 0.3109419 0.2720166 0.2719451
只有前4位数字与您输入的tol = 1e-4
相同.
Only the first 4 digits are the same, as you put tol = 1e-4
.
浮点数最多可以包含16位数字,因此建议您使用tol = 1e-16
进行可靠的收敛性测试.
A floating point number has a maximum of 16 digits, so I would suggest you use tol = 1e-16
for reliable convergence test.
这篇关于如何找到矩阵与循环收敛的时间的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!