作者empireisme (empireisme)
看板R_Language
标题Re: [问题] 回圈的计算(小复杂)
时间Sat Nov 30 23:02:04 2019
※ 引述《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://webptt.com/cn.aspx?n=bbs/R_Language/M.1575126128.A.A96.html
1F：→ a78998042a: 你两个loop中间是要赋归?12/01 10:58
2F：→ a78998042a: w = p*A*rowSums(A)12/01 10:58
3F：→ a78998042a: B1 = t(A)*p12/01 10:58
4F：→ a78998042a: w[w==0] = NA12/01 10:59
5F：→ a78998042a: w1 = 1/w12/01 10:59
6F：→ a78998042a: w1[is.na(w1)] = 012/01 10:59
7F：→ a78998042a: diag(B1%*%w1)12/01 10:59
8F：→ empireisme: 我是想要跳过分子有0的地方不要算，因为分子不为零12/01 11:12
9F：→ empireisme: 分母就不会是0 所以我才会加那个判断XD12/01 11:13
10F：→ empireisme: 甚麽是赋归?12/01 11:13
11F：→ empireisme: 我发现我新的程式码算出来的直跟我就得不一样 觉得好12/01 11:15
12F：→ empireisme: 奇怪12/01 11:16
13F：→ empireisme: 就是我只计算当alpha ij=1 我才计算 alphaij=0 不算12/01 11:35
14F：→ a78998042a: 就是你现在的写法没有在每一次loop把den归0，所以数值12/01 11:49
15F：→ a78998042a: 是迭代的，如果是迭代的，就只能用for loop12/01 11:49
我没有要迭代，只是不要算当alpha ij =0的状况，只是发现我加了这条件算出来，跟你
们不一样(你们的算法跟我旧的一样)才在想是不是我新的写错了
※ 编辑: empireisme (219.91.75.186 台湾), 12/01/2019 12:03:24
16F：→ a78998042a: 那应该是这样写 12/01 12:13
17F：→ a78998042a: for(j in 1:m){ 12/01 12:13
18F：→ a78998042a: den <- rep(0,n) 12/01 12:13
19F：→ a78998042a: densum <- rep(0,n) 12/01 12:13
20F：→ a78998042a: d<- rep(0,m) 12/01 12:13
21F：→ a78998042a: for(i in 1:n){ 12/01 12:13
22F：→ a78998042a: ... 12/01 12:13
23F：→ a78998042a: ｝ 12/01 12:13
24F：→ a78998042a: ｝ 12/01 12:13
25F：→ a78998042a: d在外面才对 12/01 12:16
※ 编辑: empireisme (219.91.75.186 台湾), 12/03/2019 18:03:30