[問題] 迴圈的計算(小複雜)

[empireisme (empireisme)]

我本身使用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://www.ptt.cc/bbs/R_Language/M.1574656447.A.C05.html
※ 編輯: empireisme (219.91.75.186 臺灣), 11/25/2019 12:39:29
→ andrew43: sum(densum)不應該和for(i)同層
→ andrew43: 你先想一下j=1;i=1時densum是什麼就能明白有問題
→ celestialgod: https://pastebin.com/TsASfKP4
→ celestialgod: 其實它就只是每個i都取代一個d[j] 但是最後一個i
→ celestialgod: 才會填一次正確值XD 是不影響結果拉...
→ celestialgod: 這個問題用sweep / colSums / rowSums就能解決了
→ celestialgod: 原po再自己看一下我上面貼的程式~ 有問題再發問
 謝謝各位回答，我是沒用過sweep,只覺得這種迴圈寫的很煩，現在沒電腦，回家再試試
看

→ andrew43: 嗯後來發現了謝謝
所以我程式應該是沒寫錯的?

因為計算量太大也不好手算qq
※ 編輯: empireisme (101.12.41.101 臺灣), 11/25/2019 14:17:10
※ 編輯: empireisme (101.12.41.101 臺灣), 11/25/2019 14:19:47
→ celestialgod: 不算寫錯 但是就多了不少操作 是可以省掉的
→ celestialgod: 這種情況建議拿 2 x 2的矩陣來驗算會好一點
→ celestialgod: sweep正如它字面上的意思就掃過去...
→ celestialgod: 自己玩玩看吧
我算過小矩陣

但就是怕有萬一qq

比如說是不是對數學公式的理解錯誤
※ 編輯: empireisme (101.12.41.101 臺灣), 11/25/2019 14:34:45
→ andrew43: 只看式子，確實多算很多。rowSums colSums先練練看
ok 看了c大的程式，看起來確實我沒有算錯但多算很多了
※ 編輯: empireisme (101.12.41.101 臺灣), 11/26/2019 10:09:21
→ a78998042a: temp = t(A)*p
→ a78998042a: den = colSums(temp)
→ a78998042a: r2 = t(temp)/den
→ a78998042a: d = colSums(r2)
→ a78998042a: 快一點，記憶體用量少一些
推 asdiy: 我的寫法是 simple=function(A,p){t=p/apply(A,1,function
→ asdiy: (x)x%*%p);apply(A,2,function(x)x%*%t} 做 10000*10000矩
→ asdiy: 陣是4秒