微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

如何在R中存储每个混合成分的观测值

如何解决如何在R中存储每个混合成分的观测值

我想手动在R中拟合混合模型。然后,我想分别存储对混合模型每个组成部分的观察。也就是说,我希望我的代码保留从每个组件淹没的观察结果。这是一个使用here中的EM算法的混合模型示例。

e_step <- function(x,mu.vector,sd.vector,alpha.vector) {


 comp1.prod <- dnorm(x,mu.vector[1],sd.vector[1]) * alpha.vector[1]
  comp2.prod <- dnorm(x,mu.vector[2],sd.vector[2]) * alpha.vector[2]
  sum.of.comps <- comp1.prod + comp2.prod
  comp1.post <- comp1.prod / sum.of.comps
  comp2.post <- comp2.prod / sum.of.comps

  sum.of.comps.ln <- log(sum.of.comps,base = exp(1))
  sum.of.comps.ln.sum <- sum(sum.of.comps.ln)

  list("loglik" = sum.of.comps.ln.sum,"posterior.df" = cbind(comp1.post,comp2.post))
}

    m_step <- function(x,posterior.df) {
  comp1.n <- sum(posterior.df[,1])
  comp2.n <- sum(posterior.df[,2])

  comp1.mu <- 1/comp1.n * sum(posterior.df[,1] * x)
  comp2.mu <- 1/comp2.n * sum(posterior.df[,2] * x)

  comp1.var <- sum(posterior.df[,1] * (x - comp1.mu)^2) * 1/comp1.n
  comp2.var <- sum(posterior.df[,2] * (x - comp2.mu)^2) * 1/comp2.n

  comp1.alpha <- comp1.n / length(x)
  comp2.alpha <- comp2.n / length(x)

  list("mu" = c(comp1.mu,comp2.mu),"var" = c(comp1.var,comp2.var),"alpha" = c(comp1.alpha,comp2.alpha))
}
for (i in 1:50) {
  if (i == 1) {
    # Initialization
    e.step <- e_step(wait,wait.summary.df[["mu"]],wait.summary.df[["std"]],wait.summary.df[["alpha"]])
    m.step <- m_step(wait,e.step[["posterior.df"]])
    cur.loglik <- e.step[["loglik"]]
    loglik.vector <- e.step[["loglik"]]
  } else {
    # Repeat E and M steps till convergence
    e.step <- e_step(wait,m.step[["mu"]],sqrt(m.step[["var"]]),m.step[["alpha"]])
    m.step <- m_step(wait,e.step[["posterior.df"]])
    loglik.vector <- c(loglik.vector,e.step[["loglik"]])

    loglik.diff <- abs((cur.loglik - e.step[["loglik"]]))
    if(loglik.diff < 1e-6) {
      break
    } else {
      cur.loglik <- e.step[["loglik"]]
    }
  }
}
loglik.vector

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。