[问题] 循环的计算(小复杂)

楼主: empireisme (empireisme)   2019-11-25 12:34:01
我本身使用R大概一两年
目前想要写出一个小function可以计算出以下的公式
https://imgur.com/supH7mE
d 是要计算出来的向量 长度是m
A 是一个矩阵维度是n*m,里面的元素都是正的
alpha是A矩阵的元素
p 是一个向量长度为m 其为机率向量 里面的元素相加等于1
以下是我的程式码,虽然检查过很多次,但不确定有没有算错
或是有没有办法用矩阵的方式去计算之,因为这样算太慢了
set.seed(1)
m=10
n=7
A <- matrix( rexp(m*n),ncol=m ,nrow=n )
p <- c( 0.05,0.1,0.3,0.05,0.05,0.05,0.07,0.03,0.2,0.1)
sum(p )
length(p )
den <- rep(0,n)
densum <- rep(0,n)
d<- rep(0,m)
for(j in 1:m){
for(i in 1:n){
den[i]<- sum(A[i,]*p)
densum[i]<- A[i,j]*p[j]/den[i]
d[j]<- sum( densum )
}
}
下面是我算出来的
[1] 0.336093550 0.855872710 2.158927311 0.233847299 0.282627585
[6] 0.253739688 0.517757929 0.312419250 0.760612933 1.288101743
作者: andrew43 (讨厌有好心推文后删文者)   2019-11-25 13:44:00
sum(densum)不应该和for(i)同层你先想一下j=1;i=1时densum是什么就能明白有问题
作者: celestialgod (天)   2019-11-25 14:09:00
https://pastebin.com/TsASfKP4其实它就只是每个i都取代一个d[j] 但是最后一个i才会填一次正确值XD 是不影响结果拉...这个问题用sweep / colSums / rowSums就能解决了原po再自己看一下我上面贴的程式~ 有问题再发问
作者: andrew43 (讨厌有好心推文后删文者)   2019-11-25 14:15:00
嗯后来发现了谢谢
作者: celestialgod (天)   2019-11-25 14:27:00
不算写错 但是就多了不少操作 是可以省掉的这种情况建议拿 2 x 2的矩阵来验算会好一点sweep正如它字面上的意思就扫过去...自己玩玩看吧
作者: andrew43 (讨厌有好心推文后删文者)   2019-11-25 14:39:00
只看式子,确实多算很多。rowSums colSums先练练看
作者: a78998042a (Benjimine)   2019-12-01 09:32:00
temp = t(A)*pden = colSums(temp)r2 = t(temp)/dend = colSums(r2)快一点,内存用量少一些
作者: asdiy (灯火阑珊)   2018-01-01 00:04:00
我的写法是 simple=function(A,p){t=p/apply(A,1,function(x)x%*%p);apply(A,2,function(x)x%*%t} 做 10000*10000矩阵是4秒

Links booklink

Contact Us: admin [ a t ] ucptt.com