Relationship of Poisson regression to determining proportions in a mixture distribution

Question

I am a particle physicist, and a very frequent task is: given data sampled from a distribution mixture $$ Z \sim \sum_i \phi_i F_i, $$ where $\phi_i$ are the mixture proportion / prior probabilities and $F_i$ are the distributions of the separate components, fit the data to measure the numbers of events from the separate components (or equivalently the weights $\phi_i$). A famous canonical example is the discovery of the Higgs boson: There are two main components here: a background whose distribution $F_b$ is exponentially falling and the Higgs signal $F_s$ whose distribution is roughly a gaussian with a mean around 125 GeV.

This class of problem is generally approached in the obvious way via likelihood maximization. Given observed values $z_k$, the log likelihood function to be maximized is $$ \ell(z_k \mid \phi) = \sum_k \log \left( \sum_i \phi_i f_i(z_k) \right) $$

I'm trying to understand the relationship of this problem to Poisson regression. My main question is: Can I perform a binned and/or unbinned mixture analysis (including the constraint $\sum_i \phi_i = 1$) using Poisson regression?

Simulation for binned case

I know that there is a relationship when the data is binned/grouped into a histogram. (Math shown at the bottom of the answer). Jumping straight to an example simulation, consider a simple case where $$ Z \sim I \mathcal{N}(0,1) + (1-I) \mathcal{N}(0,4) $$ where $I \sim B(p)$ is a latent indicator variable which is 1 with probability $p$. A sample of this model is easy to create in R:

p <- 0.7
N <- 1000
I <- rbinom(N,1,p)
Xa <- rnorm(N,0,1)
Xb <- rnorm(N,2,2)
X <- I*Xa + (1-I)*Xb

Here is the histogram of the data, hist(X):

Given this binned data, it is straightforward to estimate the number of events from the two components of the mixture using Poisson regression. The independent variables in the regression are the vectors of densities of the two mixtures at the bin centers, while the response is the vector of bin counts. The regression model uses the identity link function rather than the canonical log link, and has no intercept:

Y <- hist(X)$counts
centers <- hist(X)$mids
pa <- dnorm(centers,0,1)
pb <- dnorm(centers,2,2)
fit <- glm(Y~pa+pb-1,famil=poisson(link="identity"))

This fit does a good job estimating the mixture properties:

Call:
glm(formula = Y ~ pa + pb - 1, family = poisson(link = "identity"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6250  -0.4922   0.5906   0.7005   1.0725  

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
pa   669.61      31.53   21.24   <2e-16 ***
pb   330.86      25.63   12.91   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance:    Inf  on 13  degrees of freedom
Residual deviance: 9.6913  on 11  degrees of freedom
AIC: 78.837

Number of Fisher Scoring iterations: 4

Because there is no constraint on $\sum_i \phi_i$, the fit values seem to be counts rather than probabilities, but they appear to have the correct ratio: $669.61:330.86 = 0.67:0.33$. It could be possible to include the constraint in the glm model, replacing Y~pa+pb-1 with Y~I(pa-pb)+offset(pb)-1, but this leads to NaN errors.

I am interested in studying this analogy further, and learning what meaning the Poisson regression hat matrix has for the mixture decomposition. In particular, I'd like to know if I can extend it to an "unbinned" regression. The main conceptual difficulty I have is that in the unbinned case, the counts are all 1. The most obvious thing I can think of is

Pa <- dnorm(X,0,1)
Pb <- dnorm(X,2,2)
ufit <- glm(rep(1,N)~Pa+Pb-1,family=poisson(link="identity"))

but this has strange results:

Call:
glm(formula = rep(1, N) ~ Pa + Pb - 1, family = poisson(link = "identity"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.2657  -0.2402  -0.1172   0.1780   3.0260  

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
Pa   1.1513     0.1818   6.334  2.4e-10 ***
Pb   5.7347     0.3751  15.287  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance:    Inf  on 1000  degrees of freedom
Residual deviance: 197.91  on  998  degrees of freedom
AIC: 2201.9

Number of Fisher Scoring iterations: 6

I expect regression parameters in close to a 7:3 ratio, but they are in a nearly exactly 1:5 ratio instead.

Math details

In the binned case the histogram counts follow Poisson distributions. Letting $y_b$ be the number of events in bin $b$ and $z_b$ be the bin center, the log-likelihood function becomes $$ \ell(y_b \mid \phi) \approx \sum_b y_b \log \left( \sum_i \phi_i f_i(z_b) \right) $$ This is obviously similar to the Poisson regression log-likelihood $$ \ell_P(y_b \mid \theta) = \mathrm{constant} + \sum_b y_b \log \mu_b - \mu_b $$ where $g(\mu_b) = x_b^T \beta$ with the link function $g$. In particular, when $x_b = \left\{ f_i(z_b) \right\}$, $\beta = \phi$, and $g$ is the identity, the only difference is $$ \ell(z_b \mid \phi) - \ell_P(z_b \mid \phi) = \sum_b \mu_b = \sum_i \phi_i \left( \sum_b f_i(z_b) \right) $$ Since the $f_i(z_b)$ are pmf values (or pdf values in the unbinned case), their sum is unity, and so $$ \ell(z_b \mid \phi) - \ell_P(z_b \mid \phi) = \sum_i \phi_i $$ In the mixture decomposition problem, the weights $\phi_i$ are a priori probabilities and should sum to unity, making this difference a constant (1). I must imagine that either their is some way of incorporating this constraint in Poisson regression or that its immaterial (i.e. its presence only affects the maximum likelihood point up to a multiplicative factor $\phi' = \kappa \phi$).

Perhaps for the unbinned case, applying the constraint $\sum_i \phi_i = 1$ is more crucial. I also note that while in the binned case $\sum_b f_i(z_b) = 1$, a normalization constraint on the pmf, there is no corresponding constraint on $\sum_k f_i(z_k)$ because the values of $z_k$ are random samples -- rather than an integral over $z$ this is more like $N E[ f_i(Z)]$, the average pdf value times the number of events.