View source: R/family.mixture.R
mix2poisson  R Documentation 
Estimates the three parameters of a mixture of two Poisson distributions by maximum likelihood estimation.
mix2poisson(lphi = "logitlink", llambda = "loglink",
iphi = 0.5, il1 = NULL, il2 = NULL,
qmu = c(0.2, 0.8), nsimEIM = 100, zero = "phi")
lphi, llambda 
Link functions for the parameter 
iphi 
Initial value for 
il1, il2 
Optional initial value for 
qmu 
Vector with two values giving the probabilities relating
to the sample quantiles for obtaining initial values for

nsimEIM, zero 
See 
The probability function can be loosely written as
P(Y=y) = \phi \, Poisson(\lambda_1) +
(1\phi) \, Poisson(\lambda_2)
where \phi
is the probability an observation belongs
to the first group, and y=0,1,2,\ldots
.
The parameter \phi
satisfies 0 < \phi < 1
.
The mean of Y
is
\phi\lambda_1+(1\phi)\lambda_2
and this is returned as the fitted values.
By default, the three linear/additive predictors are
(logit(\phi), \log(\lambda_1), \log(\lambda_2))^T
.
An object of class "vglmff"
(see vglmffclass
).
The object is used by modelling functions
such as vglm
and vgam
.
This VGAM family function requires care for a successful
application.
In particular, good initial values are required because
of the presence of local solutions. Therefore running
this function with several different combinations of
arguments such as iphi
, il1
, il2
,
qmu
is highly recommended. Graphical methods such as
hist
can be used as an aid.
With grouped data (i.e., using the weights
argument)
one has to use a large value of nsimEIM
;
see the example below.
This VGAM family function is experimental and should be used with care.
The response must be integervalued since
dpois
is invoked.
Fitting this model successfully to data can be difficult
due to local solutions and illconditioned data. It pays to
fit the model several times with different initial values,
and check that the best fit looks reasonable. Plotting
the results is recommended. This function works better as
\lambda_1
and \lambda_2
become
more different. The default control argument trace =
TRUE
is to encourage monitoring convergence.
T. W. Yee
rpois
,
poissonff
,
mix2normal
.
## Not run: # Example 1: simulated data
nn < 1000
mu1 < exp(2.5) # Also known as lambda1
mu2 < exp(3)
(phi < logitlink(0.5, inverse = TRUE))
mdata < data.frame(y = rpois(nn, ifelse(runif(nn) < phi, mu1, mu2)))
mfit < vglm(y ~ 1, mix2poisson, data = mdata)
coef(mfit, matrix = TRUE)
# Compare the results with the truth
round(rbind('Estimated' = Coef(mfit), 'Truth' = c(phi, mu1, mu2)), 2)
ty < with(mdata, table(y))
plot(names(ty), ty, type = "h", main = "Orange=estimate, blue=truth",
ylab = "Frequency", xlab = "y")
abline(v = Coef(mfit)[1], lty = 2, col = "orange", lwd = 2)
abline(v = c(mu1, mu2), lty = 2, col = "blue", lwd = 2)
# Example 2: London Times data (Lange, 1997, p.31)
ltdata1 < data.frame(deaths = 0:9,
freq = c(162,267,271, 185,111,61,27,8,3,1))
ltdata2 < data.frame(y = with(ltdata1, rep(deaths, freq)))
# Usually this does not work well unless nsimEIM is large
Mfit < vglm(deaths ~ 1, weight = freq, data = ltdata1,
mix2poisson(iphi=0.3, il1=1, il2=2.5, nsimEIM=5000))
# This works better in general
Mfit = vglm(y ~ 1, mix2poisson(iphi=0.3, il1=1, il2=2.5), ltdata2)
coef(Mfit, matrix = TRUE)
Coef(Mfit)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.