Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Kirin_Guess
    started a topic Sum of an infinite series

    Sum of an infinite series

    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.
    Last edited by Kirin_Guess; 06 Feb 2019, 18:43.

  • Jean-Claude Arbaut
    replied
    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
    Last edited by Jean-Claude Arbaut; 26 Apr 2019, 04:01.

    Leave a comment:


  • Jean-Claude Arbaut
    replied
    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.

    Code:
    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

    Leave a comment:


  • Jeronimo Muniz
    replied
    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.

    Leave a comment:

Working...
X