Hello there. I would like to know and learn how Stata staff usually does when programming a statistical method that involves the sum of an infinite series. Suppose that each subject i in a sample of size n has two random variables: a latent variable Yi follows a Poisson distribution with mean λi, and an observed variable Xi follows a Bernoulli distribution with probability πi that is certain function of the first variable, i.e. πi = f(Yi). Does Stata (Mata in particular) provide any function(s) for approximating the following summation?

\[ \sum_{y=0}^{\infty} f(Y_i)^{X_i} [1-f(Y_i)]^{1-X_i} \biggl( \frac{\lambda_i^{Y_i}}{Y_i!}e^{-\lambda_i} \biggl) \]

The only way I know for the moment is to use the while loop and do the summation from y=0 and stop the loop when the change between iterations is smaller than a certain value, e.g. 1e-8. However, this way is inefficient, because subjects have different λi and Xi, and the numbers of loop to achieve the stop criterion vary over subjects. Therefore I am looking for more clever methods.

I know that, for certain forms of continous variables, Stata performs Gauss–Hermite quadrature to approximate the values of integrals, and the undocumented Mata function _gauss_hermite_nodes() can help us to do so. I wonder whether there are Stata/Mata functions that can help out in approximating the sum of an infinite series for discrete variabels.

Thanks.