0

我有一个简单的泊松模型，具有未观察到的正态分布的异质性 eta。这可以通过将 eta 从对数似然中积分来解决。但是，当我尝试这样做时，它会产生错误的系数，我不明白为什么。代码如下所示。它应该相应地返回接近 0.018 和 0.013 的 beta 和 sigma 系数，但事实并非如此。

library("maxLik")
library("pracma")
library("rmutil")

n = 3000
x = rnorm(n, mean=0, sd=2)
beta1 = 0.018
sigma = 0.013
eta = rnorm(n, mean=0, sd= sigma)
mu = exp(x*beta1 + eta)
y = rpois(n,mu)


poisson <- function(y,mu) {
  a = mu^y*exp(-mu)/factorial(y)
  return(a)
}




loglik <- function(param) {
  beta = param[1]
  sigma = param[2]
  f <- function(eta) dnorm(eta, sd = sigma)*poisson(y, exp((x*beta + eta)))
  ll <- int(f, a = -100, b = 100) 
  ll[ll==0] <- 0.00000000001
  logll = log(ll)
  return((sum(logll)))
}
bet = runif(16, min=-1, max = 0)/10
ml <- maxLik(loglik, start = c(beta= 0.1, sigma = 0.1), print.level = 2)
4

0 回答 0