从网上搜了半天,好像只有 %*%可以做矩阵乘法。
如果我想算一个矩阵的4,5,6次方,我应该怎么操作?试了试直接^4好像不行(记得matlab里面好像是这样的),只是把矩阵里的每个数分别算了个四次方。
问题白痴,还请别见怪。或许R不应该用来计算这种东西,但是我觉得应该能算吧?懒的再装matlab了。
问个白痴问题——矩阵的N次方怎么算?
最笨的方法,呵呵:
<br />
X%*%X ##二次方<br />
X%*%X*%X ##三次方<br />
X%*%X%*%X%*%X ##四次方<br />
</p>
比较赞、值得收藏的方法:http://cos.name/cn/topic/5619
根据 http://bm2.genes.nig.ac.jp/RGM2/R_current/library/Biodem/man/mtx.exp.html
<br />
> library(Biodem)<br />
> test<-matrix(c(1:16), 4,4)<br />
> pow.test<-mtx.exp(test,10)<br />
> pow.test<br />
[,1] [,2] [,3] [,4]<br />
[1,] 2.568153e+14 5.933637e+14 9.299121e+14 1.266461e+15<br />
[2,] 2.908200e+14 6.719305e+14 1.053041e+15 1.434152e+15<br />
[3,] 3.248247e+14 7.504974e+14 1.176170e+15 1.601843e+15<br />
[4,] 3.588295e+14 8.290642e+14 1.299299e+15 1.769534e+15</p>
<p>
</p>……我想了半天,搜了立方、N次方、power等等,就没想到“乘方”,汗。
yp版主的解决办法真赞,这个函数挺好用的,嘿嘿。对我这种懒人来说,拷贝一个程序不如下个包来得快。
顺便说一下在看到这个帖子前我的解决办法……别见笑。
<br />
b <- a %*% a<br />
b <- b%*%b<br />
然后我就不停的按上箭头和回车重复……汗。其实主要是想看一个转移矩阵的长期的稳定值……后来发现干脆用手算解方程,再汗。
</p>
回复 第4楼 的 Cloudly:
其实可以用递归的方法构造矩阵乘方的函数
<br />
mpower <- function(mat,n){<br />
if(n == 1)<br />
return(mat)<br />
else<br />
return(mat %*% mpower(mat,n-1))<br />
}<br />
</p>
但是我觉得如果 n很大的话 一来速度会变慢 二来误差也会积累下去
刚才下载了Biodem的包 看到mtx.exp的代码 原来也这么简单 循环比递归快
<br />
mtx.exp <- function (X, n)<br />
{<br />
if (n != round(n)) {<br />
n <- round(n)<br />
warning("rounding exponent `n' to", n)<br />
}<br />
phi <- diag(nrow = nrow(X))<br />
pot <- X<br />
while (n > 0) {<br />
if (n%%2)<br />
phi <- phi %*% pot<br />
n <- n%/%2<br />
pot <- pot %*% pot<br />
}<br />
return(phi)<br />
}</p>
<p>
</p>
PS: 转移矩阵是是对称矩阵 它的乘方是有规律可寻的
该矩阵可以分解成 A=LDL' D是个对角阵
A^n = L(D^n)L' 而对角阵的n次方也还是对角阵 对应元素做乘方就行
回复 第5楼 的 ypchen:……偏偏俺遇到的不是个对称阵……可能没说明,是马尔可夫转移矩阵,各种状态之间转移的……所以不一定会是对角阵吧?(状态A->B的概率一定等于B->A的概率?)
其实马尔可夫转移矩阵的话,手动解方程就能求长期稳定状态了,汗。我开始只是不想解方程,就想算算N次方,后来还是乖乖的解方程去了(其实貌似也可以用R解)……
至于递归和循环,指数和线性的关系么?