4
$\begingroup$

I am trying to fit a model to some simulated data. The idea is to use ML-Estimation. However, I am totally lost.

I have a dependent variable which is a sum of two (unknown) variables. The first part is Poisson-Distributed and the second part is Negative Binomial Distributed. In the sum this leads to a Delaporte model.

I already implemented parts of the model based on a log likelihood-function. However, I am missing the weighting matrix.

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)
out

Where $???$ is needed! $$ ??? = \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}

$\endgroup$
0
1
$\begingroup$

Have you considered using the Delaporte package for R? Try using optim or nloptr to minimize the NLL as returned by ddelap in that package. It may take some time as the Delaporte has no closed form and requires summation.

For example:

set.seed(6541321)    
library(nloptr)
library(Delaporte)
POIS <- rpois(5000, lambda = 4)
summary(POIS)
>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    3.00    4.00    4.03    5.00   13.00 
NB <- rnbinom(5000, size = 3, mu = 15) ## Implies alpha = 3 and beta = 5
summary(NB)
>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    8.00   13.00   15.29   20.00   81.00 
Del_Test <- POIS + NB
summary(Del_Test)
>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   12.00   17.00   19.32   25.00   86.00 
Inits <- MoMdelap(Del_Test)
Inits ## Method of Moments estimate (if reasonable)
>[1] 2.777264 5.337955 4.497891
Del_NLL <- function(pars, data) {
  return(-sum(ddelap(data, pars[[1]], pars[[2]], pars[[3]], log = TRUE)))
}
Del_Fit <- nloptr(x0 = Inits, eval_f = Del_NLL, data = Del_Test, lb = c(0, 0, 0), opts=list(algorithm = "NLOPT_LN_SBPLX", maxeval = 10000))
Del_Fit
>Call:
nloptr(x0 = Inits, eval_f = Del_NLL, lb = c(0, 0, 0), opts = list(algorithm = "NLOPT_LN_SBPLX",     maxeval = 10000), data = Del_Test)


Minimization using NLopt version 2.4.0 

NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached. )

Number of Iterations....: 2250 
Termination conditions:  maxeval: 10000 
Number of inequality constraints:  0 
Number of equality constraints:    0 
Optimal value of objective function:  18061.4279834719 
Optimal value of controls: 3.026119 5.1027 3.882561
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.