First of all, I am rather new to statistics, so go easy on me.
I am aware that the negative binomial distribution can be thought to arise as a result of letting the $\lambda$ parameter in a Poisson distribution vary like the Gamma distribution. I wanted to know if there was a similar result for a Poisson-Gaussian mix, where $\lambda$ was a random variable and was distributed as a Gaussian. Obviously there is the issue of Gaussian distributions allowing for negative values, but if we restrict it from zero to infinity and normalize accordingly I thought it should be possible.
In this case, the $\lambda$ parameter would have some average value and a symmetric spread around that value.
From what I understand, the mixture is represented by following integral.
$$\frac{\int_{0}^{\infty}\frac{\lambda^k}{k!}e^{-\lambda} e^{-\frac{(\lambda-\mu)^2}{2 \sigma^2}}d\lambda}{\int_{0}^{\infty}e^{-\frac{(\lambda-\mu)^2}{2 \sigma^2}}d\lambda}$$
If a simple or well-known form does not exist, how would you all recommend implementing a fit for this in a program like Matlab? Computing it from this definition and fitting to it seems rather difficult. I was also looking into a mixing with a log-normal distribution (since $\lambda$ is positive), but I thought that would be even more difficult to deal with.