Fitting Mixture of Poissons Without GLM

Question

I am looking for a way to fit a mixture of 2 univariate Poisson distributions, but the R packages I checked, mixtools and flexmix, assume a GLM model. They require a "formula", etc.. How can I fit a non-GLM model to my data?

In Hidden Markov Models for Time Series: An Introduction Using R by MacDonald et al. one exercise suggests to fit Poisson distributions using nlm of optim; I guess that is possible too, the error function could be tied to the likelihood of the data, the routine would optimize on that.

Note: Per @Tim's comment, I used

library("flexmix")
y <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,
       1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,2,1,1,1,
       1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,
       0,1,0,0,0,2,1,0,0,0,1,1,0,2,3,3,1,1,2,1,1,1,
       1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1)
df <- data.frame(y)
res <- flexmix(y ~ 1, data = df, k = 2, model = FLXMRglm(family = "poisson"))
print (summary(res))
print (posterior(res))

This seemed to work; my goal was fitting two Poisson distributions to coal mining accident data to find a switch point between two different regimes.

Thanks,

What is the problem with GLM and formula for you? You can always use `variable ~ 1` formula for intercept-only (i.e. mean only) model... — Tim, Nov 13 '15 at 09:01
It's not a problem per se - I dont need it. I took your advice, and it seemed to work. I will update the question. — BBSysDyn, Nov 13 '15 at 09:12
Glad it helped, I posted a more detailed answer with some references FYI. — Tim, Nov 13 '15 at 10:06
Another quick question: the model above returns me 1.0101464 and -0.580075 - my per-cluster Poisson lambda's? One value is negative though, aren't Poisson lambdas supposed to be all positive? — BBSysDyn, Nov 13 '15 at 13:27
Well... Poisson's $\lambda$ are only non-negative integers, but you are using GLM with Poisson link function - that is not the same. Take a look at literature related to GLM and multiple threads here on GLM and link functions like http://stats.stackexchange.com/questions/20523/difference-between-logit-and-probit-models#30909 to understent them better. — Tim, Nov 13 '15 at 13:51

score 3 · Accepted Answer · edited Apr 13 '17 at 12:44

You can use any GLM model (see more about GLM in here) in univariate case, if the general case is

$$ Y = \beta_0 + \beta_1 X + \varepsilon $$

then you can use intercept-only model

$$ Y = \beta_0 + \varepsilon $$

(or in R formula Y ~ 1). Such model simply estimates the mean, e.g.

> mean(mtcars$mpg)
[1] 20.09062
> lm(mpg ~ 1, mtcars)

Call:
lm(formula = mpg ~ 1, data = mtcars)

Coefficients:
(Intercept)  
      20.09

This can be extended into other cases, like in the one you described. In fact, if you look deeper into flexmix documentation (see also both JSS papers here and here and multiple papers by Bettina Grün and Friedrich Leisch that are available online), you'll see that multiple examples deal with such intercept-only formulas where the FLX... part is used for more advanced features of the model.

score 0 · Answer 2 · answered Nov 20 '15 at 11:00

Two code examples I found, one based on R nlm, the other Matlab's mle.

R

xf <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,2,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,3,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1)

udist <- function(n) rep(1/n, n) 

#natural to working
n2wp <- function(p) {
    m <- length(p)
    log(p[2:m]/(1 - sum(p[2:m])))
}

#working to natural
w2np <- function(lp) {
    rv <- exp(lp)/(1 + sum(exp(lp)))    
    c(1 - sum(rv), rv)
}

#optimisation function
of <- function(pv, m, x) {
    #convert working parameters to natural paramters
    pr <- exp(pv[1:m])
    probs <- w2np(pv[(m+1):(2*m - 1)])
    #calculate -ve log likelihood
    -sum(log(outer(x, pr, dpois) %*% probs))
}

#initial estimates and probabilities for 2, 3 and 4 distributions
#the lambda values I just guess, and use an uniform distribution
#for the initial mixing distribution.
pv <- c(log(c(1, 2)), n2wp(udist(2)))

#number of distributions to fit
m <- 2

#fit using nlm
fv <-nlm(of, pv, m, xf, print.level=0) 
rv <- fv$est

#lambda estimates
exp(rv[1:m])
#mixing distribution
w2np(rv[(m+1):(2*m-1)])

Matlab

x = [poissrnd(1,100,1); poissrnd(10,200,1)];
f = @(x,p,logmu1,logmu2) p*poisspdf(x,exp(logmu1)) + (1-p)*poisspdf(x,exp(logmu2));
mle(x,'pdf',f,'start',[.5 2 4])
exp(ans(2:3))

Fitting Mixture of Poissons Without GLM

2 Answers2