Fitting Mixture of Poissons Without GLM 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,
 A: 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.
A: 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))

