作者empireisme (empireisme)
看板R_Language
标题[问题] 回圈的计算(小复杂)
时间Mon Nov 25 12:34:01 2019
我本身使用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
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 219.91.75.186 (台湾)
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/R_Language/M.1574656447.A.C05.html
※ 编辑: empireisme (219.91.75.186 台湾), 11/25/2019 12:39:29
1F：→ andrew43: sum(densum)不应该和for(i)同层11/25 13:44
2F：→ andrew43: 你先想一下j=1;i=1时densum是什麽就能明白有问题11/25 13:44
4F：→ celestialgod: 其实它就只是每个i都取代一个d[j] 但是最後一个i11/25 14:11
5F：→ celestialgod: 才会填一次正确值XD 是不影响结果拉...11/25 14:11
6F：→ celestialgod: 这个问题用sweep / colSums / rowSums就能解决了11/25 14:12
7F：→ celestialgod: 原po再自己看一下我上面贴的程式~ 有问题再发问11/25 14:12
谢谢各位回答，我是没用过sweep,只觉得这种回圈写的很烦，现在没电脑，回家再试试
看
8F：→ andrew43: 嗯後来发现了谢谢11/25 14:15
所以我程式应该是没写错的?
因为计算量太大也不好手算qq
※ 编辑: empireisme (101.12.41.101 台湾), 11/25/2019 14:17:10
※ 编辑: empireisme (101.12.41.101 台湾), 11/25/2019 14:19:47
9F：→ celestialgod: 不算写错 但是就多了不少操作 是可以省掉的11/25 14:27
10F：→ celestialgod: 这种情况建议拿 2 x 2的矩阵来验算会好一点11/25 14:28
11F：→ celestialgod: sweep正如它字面上的意思就扫过去... 11/25 14:29
12F：→ celestialgod: 自己玩玩看吧11/25 14:29
我算过小矩阵
但就是怕有万一qq
比如说是不是对数学公式的理解错误
※ 编辑: empireisme (101.12.41.101 台湾), 11/25/2019 14:34:45
13F：→ andrew43: 只看式子，确实多算很多。rowSums colSums先练练看11/25 14:39
ok 看了c大的程式，看起来确实我没有算错但多算很多了
※ 编辑: empireisme (101.12.41.101 台湾), 11/26/2019 10:09:21
14F：→ a78998042a: temp = t(A)*p 12/01 09:32
15F：→ a78998042a: den = colSums(temp) 12/01 09:32
16F：→ a78998042a: r2 = t(temp)/den 12/01 09:32
17F：→ a78998042a: d = colSums(r2) 12/01 09:33
18F：→ a78998042a: 快一点，记忆体用量少一些 12/01 09:33
19F：推 asdiy: 我的写法是 simple=function(A,p){t=p/apply(A,1,function 01/01 00:04
20F：→ asdiy: (x)x%*%p);apply(A,2,function(x)x%*%t} 做 10000*10000矩 01/01 00:04
21F：→ asdiy: 阵是4秒 01/01 00:04