Sum of an infinite series

[Kirin_Guess]

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 Y_i follows a Poisson distribution with mean λ_i, and an observed variable X_i follows a Bernoulli distribution with probability π_i that is certain function of the first variable, i.e. π_i = f(Y_i). 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 X_i, 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.

[Jeronimo Muniz]

A similar challenge would be to solve the following equation: A=exp(nM), in which A and M art square matrices of the same order. In this equation, exp(nM) defines an infinite power series, and A is known. How should one write a Mata program to solve for M? Help is deeply appreciated.

[Jean-Claude Arbaut]

Jeronimo Muniz

This should be another question. The matrix exponential is unique, but the matrix logarithm is not. However, it's fairly easy to find a solution, using the eigendecomposition.

mata
a=rnormal(4,4,0,1)

/* compute the exponential */
eigensystem(a,p=.,d=.)
b=p*diag(exp(d))*qrinv(p)

/* compute the logarithm */
eigensystem(b,p=.,d=.)
c=p*diag(log(d))*qrinv(p)

norm(a-c)
end

There are other ways. For instance, when the matrix is not diagonalizable, you can use the Lagrange-Hermite interpolation polynomial on the matrix spectrum (see this).
See also "Functions of Matrices" by Nick Higham.

HTH

Jean-Claude Arbaut

[Jean-Claude Arbaut]

Kirin_Guess

This doesn't answer the question, but note that _gauss_hermite_nodes is actually documented. However, _gauss_laguerre_nodes and _gauss_legendre_nodes are not.

Jean-Claude Arbaut