Estoy intentando ajustar un modelo a unos datos simulados. La idea es utilizar la estimación ML. Sin embargo, estoy totalmente perdido.
Tengo una variable dependiente que es una suma de dos variables (desconocidas). La primera parte tiene una distribución de Poisson y la segunda tiene una distribución binomial negativa. En la suma esto lleva a un Modelo Delaporte .
Ya he implementado partes del modelo basadas en una función de probabilidad logarítmica. Sin embargo, me falta la matriz de ponderación.
set.seed(1234)
x1 <- runif(5000, min = 0, max = 5)
x2 <- runif(5000, min = 0, max = 5)
# Model 1
## beta0=0.5, beta1=0.3, beta2=0.4, eta=0.5, delta=2
alpha1 <- rgamma(5000, shape = 0.5, rate = 0.5)
lambda1 <- alpha1 * exp(0.5 + 0.3 * x1 + 0.4 * x2 )
Q1 <- matrix(1:5000, ncol=1)
for (i in 1:5000){
Q1[i] <- rpois(1, lambda=lambda1)
}
Q1 <- apply(Q1, 1,as.numeric) # converts matrix object to numeric vector
u1 <- rpois(5000, lambda=0.5)
y1 <- Q1 + u1
# ML Estimation
n <- length(x1)
nlogL <- function(par) {
n <- length(x1)
beta <- par[1:3]
eta <- par[4]
delta <- par[5]
-(-n*eta*log(eta)*n*mean(y1)+sum(delta*log(delta/(delta+lambda)))+sum(log(???)))
}
par0 <- as.vector(c(0.2, 0.1, 0.2, 0.2, 1.3))
out<-nlm(nlogL,p=par0, hessian = TRUE)
outDonde $???$ ¡es necesario! $$ ??? = \sum_{q_i = 1}^{y_i} w(q_i) $$ \begin{align} w &= w(q_i) \\ &= w(q_i| \beta, \eta, \delta) \\ &= \frac{\Gamma(q_i + \delta)}{\Gamma(\delta)\Gamma(Q_i+1} \bigg[\frac{\lambda_i}{\eta(\lambda_i+\delta)}\bigg]^q_i \frac{1}{\Gamma(y_1-q_i+1)} \end{align}
