Normal-poisson composite density

Modeling
1

Hi

I am trying to fit a model to data generated as counts but measured with normal error - e.g. a particle detector which detects the particles using an analogue detector:

count ~poisson(lambda)
measured_signal ~ normal(count,sd)

This runs into problems due to the parameter counts being an integer. I understand a way around this would be to analytically sum over parameter ‘counts’ to find the density of the compound distribution, but cant work out how to do this, nor could mathmatica. Any thoughts for a way round this?

simulate data (in r):
#simulate data
lambda = 2
sigma = .5
counts = rpois(1e3,lambda)
measure = rnorm(1e3,counts,sigma)
hist(measure)

Stan model (wont work as counts is int)

model {
data {
  int N;
  real measure[N];
}
parameters {
  real<lower=0> lambda;
  int<lower=0> counts[N];
  real<lower=0> sigma;
}
model {
counts ~ poisson(lambda);
measure ~ normal(counts,sigma);
}
}
2

Depends on how large count and sigma are. If sigma<1 then you can usually recover count just by rounding the measured value to the nearest integer. And if count is very large then poisson is approximately normal so that’s easy too. But I assume you’re not dealing with these easy scenarios.

Numerical summation is still possible if there are not too many terms. I’d try something like this

data {
  int N;
  real measure[N];
}
parameters {
  real<lower=0> lambda;
  int<lower=0> counts[N];
  real<lower=0> sigma;
}
model {
  for (i in 1:N) {
    if (measure[i] < 200) {
      // count is probably small, direct summation is feasible
      real ll[301];
      for (count in 0:300) // count > 300 has negligible probability
        ll[1+count] = poisson_lpmf(count| lambda) + normal_lpdf(measure[i]| count, sigma);
      target += log_sum_exp(ll);
    } else {
      // count is large, normal approximation may be adequate
      target += normal_lpdf(measure[i]| lambda, sqrt(lambda + square(sigma)));
    }
  }
}
2 Likes
3

You will need to separate out the measured zeroes in the above solution. You might also check out the Tweedie distribution which is a compound Poisson gamma, which could be close enough to what you want.

4

Another trick would be to approximate Poisson by gamma distribution by matching the first and second moments (then the parameters would be shape=lambda and invscale=1)