PTT數位生活區 / R_Language

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

看板R_Language作者 (empireisme)時間4年前 (2019/11/30 23:02), 4年前編輯推噓0(0推 0噓 25→)
留言25則, 2人參與, 4年前最新討論串2/2 (看更多)
※ 引述《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
→
12/01 10:58, 4年前 , 1F
你兩個loop中間是要賦歸?
12/01 10:58, 1F
→
12/01 10:58, 4年前 , 2F
w = p*A*rowSums(A)
12/01 10:58, 2F
→
12/01 10:58, 4年前 , 3F
B1 = t(A)*p
12/01 10:58, 3F
→
12/01 10:59, 4年前 , 4F
w[w==0] = NA
12/01 10:59, 4F
→
12/01 10:59, 4年前 , 5F
w1 = 1/w
12/01 10:59, 5F
→
12/01 10:59, 4年前 , 6F
w1[is.na(w1)] = 0
12/01 10:59, 6F
→
12/01 10:59, 4年前 , 7F
diag(B1%*%w1)
12/01 10:59, 7F
→
12/01 11:12, 4年前 , 8F
我是想要跳過分子有0的地方不要算，因為分子不為零
12/01 11:12, 8F
→
12/01 11:13, 4年前 , 9F
分母就不會是0 所以我才會加那個判斷XD
12/01 11:13, 9F
→
12/01 11:13, 4年前 , 10F
甚麼是賦歸?
12/01 11:13, 10F
→
12/01 11:15, 4年前 , 11F
我發現我新的程式碼算出來的直跟我就得不一樣 覺得好
12/01 11:15, 11F
→
12/01 11:16, 4年前 , 12F
奇怪
12/01 11:16, 12F
→
12/01 11:35, 4年前 , 13F
就是我只計算當alpha ij=1 我才計算 alphaij=0 不算
12/01 11:35, 13F
→
12/01 11:49, 4年前 , 14F
就是你現在的寫法沒有在每一次loop把den歸0，所以數值
12/01 11:49, 14F
→
12/01 11:49, 4年前 , 15F
是迭代的，如果是迭代的，就只能用for loop
12/01 11:49, 15F
我沒有要迭代，只是不要算當alpha ij =0的狀況，只是發現我加了這條件算出來，跟你 們不一樣(你們的算法跟我舊的一樣)才在想是不是我新的寫錯了 ※ 編輯: empireisme (219.91.75.186 臺灣), 12/01/2019 12:03:24
→
12/01 12:13, 4年前 , 16F
那應該是這樣寫
12/01 12:13, 16F
→
12/01 12:13, 4年前 , 17F
for(j in 1:m){
12/01 12:13, 17F
→
12/01 12:13, 4年前 , 18F
den <- rep(0,n)
12/01 12:13, 18F
→
12/01 12:13, 4年前 , 19F
densum <- rep(0,n)
12/01 12:13, 19F
→
12/01 12:13, 4年前 , 20F
d<- rep(0,m)
12/01 12:13, 20F
→
12/01 12:13, 4年前 , 21F
for(i in 1:n){
12/01 12:13, 21F
→
12/01 12:13, 4年前 , 22F
...
12/01 12:13, 22F
→
12/01 12:13, 4年前 , 23F
｝
12/01 12:13, 23F
→
12/01 12:13, 4年前 , 24F
｝
12/01 12:13, 24F
→
12/01 12:16, 4年前 , 25F
d在外面才對
12/01 12:16, 25F
※ 編輯: empireisme (219.91.75.186 臺灣), 12/03/2019 18:03:30
‣ 更多 empireisme 的文章
文章代碼(AID): #1TueHmgM (R_Language)
討論串 (同標題文章)
本文引述了以下文章的的內容：
完整討論串 (本文為第 2 之 2 篇)：
排序： 最舊先 | 最新先 | 留言數
在新視窗開啟完整討論串 (共2篇)
R_Language 近期熱門文章
更多 近期熱門文章 >>
PTT數位生活區 即時熱門文章
更多 即時熱門文章 >>
‣ 更多 empireisme 的文章
文章代碼(AID): #1TueHmgM (R_Language)