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

[empireisme (empireisme)]

※ 引述《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

https://pastebin.com/wBRDricf

我需要計算的A矩陣有0在裡面，這是我後來的想法，關鍵點在於

我新增了 一行A(i,j)>0的判斷

不知道如果我的A矩陣 有很多0的狀況下 是否就只能用迴圈來寫了

或是有更快的方法 謝謝

我新增的寫法
確認這方法可以用 所以補上

expected_death3 <- function(A,s){
#A <- datatest_A
#s <- s_test
numerator<- sweep(A, 2, s, "*")
#numerator
denominator <- matrix(rep(A%*%s,length(s)),ncol=length(s))
denominator
bigA <- numerator/denominator
#bigA
bigA[is.nan(bigA)] <- 0 #to skip zero
#bigA
d <- colSums(bigA)
return(d)
}
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 219.91.75.186 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1575126128.A.A96.html
→ a78998042a: 你兩個loop中間是要賦歸?
→ a78998042a: w = p*A*rowSums(A)
→ a78998042a: B1 = t(A)*p
→ a78998042a: w[w==0] = NA
→ a78998042a: w1 = 1/w
→ a78998042a: w1[is.na(w1)] = 0
→ a78998042a: diag(B1%*%w1)
→ empireisme: 我是想要跳過分子有0的地方不要算，因為分子不為零
→ empireisme: 分母就不會是0 所以我才會加那個判斷XD
→ empireisme: 甚麼是賦歸?
→ empireisme: 我發現我新的程式碼算出來的直跟我就得不一樣 覺得好
→ empireisme: 奇怪
→ empireisme: 就是我只計算當alpha ij=1 我才計算 alphaij=0 不算
→ a78998042a: 就是你現在的寫法沒有在每一次loop把den歸0，所以數值
→ a78998042a: 是迭代的，如果是迭代的，就只能用for loop
我沒有要迭代，只是不要算當alpha ij =0的狀況，只是發現我加了這條件算出來，跟你
們不一樣(你們的算法跟我舊的一樣)才在想是不是我新的寫錯了
※ 編輯: empireisme (219.91.75.186 臺灣), 12/01/2019 12:03:24
→ a78998042a: 那應該是這樣寫
→ a78998042a: for(j in 1:m){
→ a78998042a:     den <- rep(0,n)
→ a78998042a:     densum <- rep(0,n)
→ a78998042a:     d<- rep(0,m)
→ a78998042a: for(i in 1:n){
→ a78998042a: ...
→ a78998042a: ｝
→ a78998042a: ｝
→ a78998042a: d在外面才對
※ 編輯: empireisme (219.91.75.186 臺灣), 12/03/2019 18:03:30