# Fitting Mixture Distribution of Poisson and Negative Binomial (Delaporte Distribution)

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}

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