R数值计算多重求和编程

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

        附件的图片是要解决的问题，目标是对最下行的极大似然函数进行估计，其中的u，k，w是三个未知参数，需要估计出来。t是时间，T是最后一次事件发生的时间，t下标n是第n次事件发生的时间，t下标k是第k次事件发生的时间。左侧的对数求和项包括两个累加求和，现在尝试用了nlminb和optim两个函数，但得到的结果都不理想。我把这个式子直接解读为：

  c<- sum(log(u+k*w*sum(exp(-1*w*m^2/t))))-u*T-k*sum(t*exp(w*m*(T-t)/t-1))

        红色的sum是想实现对每个t，求t(1)-t到t(n-1)-t的指数和；绿色的sum是对每个取值t求对数和；

        其中的m是各项t的tk-ti。

        这个式子是错误的，正确的式子似乎应该包括至少两个循环。并且如何在参数估计的公式里加入循环，我也不懂，网上的教程都是针对简单案例的，对解决我的问题帮助不大。我对R编程还不是很了解，希望能在大家的帮助下解决这个问题。谢谢！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：数值计算 nlminb Optim 极大似然 似然函数 图片

# 创建一个向量T, T中包括m个时间，每个时间就是你帖子中的时间k
T <- c(all the values in vector t, the number of elements in vector t should be m, and m is the loop length in k)
inner.sum <- 0   # 初始化inner.sum，用来对第二个for loop里的内容求和，也就是你贴子中红色的sum
outer.sum <- 0  # 初始化outer.sum，用来对第一个for loop里的内容求和，也就是你帖子中绿色的sum
for (t in 1:n) {
  for (k in 1:m) {
    if (T[k] >= t) {
      next
    } else {
      inner.sum <- inner.sum + w * exp(-w * (t - T[k]))  #对红色sum累加求和
    }
   
  }
  k0_inner.sum <- k0 * inner.sum  
  log_inner.sum <- ln(u + k0_inner.sum)  
  outer.sum <- out.sum + log_inner.sum  #对绿色sum累加求和
}

print (outer.sum)

这是个大概的code，不包括你公式中的u*T-k*sum(t*exp(w*m*(T-t)/t-1))。

zhgzhckc 发表于 2015-3-25 05:14
# 创建一个向量T, T中包括m个时间，每个时间就是你帖子中的时间k
T

非常感谢！我先试一下。

zhgzhckc 发表于 2015-3-25 05:14
# 创建一个向量T, T中包括m个时间，每个时间就是你帖子中的时间k
T

您好，按照您的思路，我调整了代码，已经可以用来解决问题了，非常感谢。下面是我用来做极大似然估计的代码，在运行小样本的时候，耗时还可以接受，n取200左右，需要5秒左右的响应时间可以出结果。

现在有新的问题，当n或者m较大，比如有30000到50000个时，有没有办法能够加速运算？我加入数据算了5个小时还没出现结果。。。

只能再向大神求教了！

LL<- function(x)
{
  u<- x[1]
  w<- x[2]
  k<- x[3]
  t<- with(Data, count)
  inner.sum<-0
  outer.sum<-0
  extra.sum<-0
  for (count in 1:33489) {
    for (b in 1:33489) {
      if (t >= count) {
        next
      } else {
        inner.sum <- inner.sum + w * exp(-w * (count - t))
      }
        extra.sum <- extra.sum + exp(w * t)-exp(-w * (4355 - t))      
    }
    k0_inner.sum <- k * inner.sum  
    log_inner.sum <- log(u + k0_inner.sum)  
    outer.sum <- outer.sum + log_inner.sum
    k0_extra.sum <- k * extra.sum
    c<- outer.sum-u*4355-k0_extra.sum   
  }
  return(-c)
}
hist(count,freq=F)
lines(density(count))
Data.res<- nlminb(c(0.001,0.001,0.001),LL,
                  lower=c(0.0001,0.0001,0.0001),
                  upper=c(0.9999,0.9999,0.9999))
Data.res$par