Announcement

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

  • 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.

  • #2
    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.

    Comment


    • #3
      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

      Comment


      • #4
        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.

        Comment

        Working...
        X