hgm R

Software Packages for Holonomic Gradient Method. The numerical evaluation of the normalizing constant for a given statistical distribution is a fundamental problem in statistics. For example, the normalizing constant of the Gaussian distribution is expressed in terms of a rational expression of a parameter of the distribution named as the standard deviation. However, normalizing constants of many interesting stasistical distributions do not have such closed expressions. The holonomic gradient method, HGM in short, is a general method to evaluate normalizing constant numerically for several parameters in the framework of Zeilberger’s holonomic systems approach. In fact, broad classes of normalizing constants are holonomic functions with respect to parameters. Then, such normalizing constants satisfy holonomic systems of linear partial differential equations. The HGM consists of three steps for a given normalizing constant. (1) Find a holonomic system satisfied by the normalizing constant. We may use computational or theoretical methods to find it. Groebner basis and related methods are used. (2) Find an initial value vector for the holonomic system. This is equivalent to evaluating the normalizing constant and its derivatives at a point. This step is usually performed by a series expansion. (3) Solve the holonomic system numerically. We use several methods in numerical analysis such as the Runge-Kutta method of solving ordinary differential equations and efficient solvers of systems of linear equations. The HGM was proposed in 2011 by a group of people inclusing us and has given several new results. For example, the orthant probability is the normalizing constant of the multivariate normal distribution restricted to the first orthant. The HGM can evaluate it in a high accuracy up to the 20 dimensional case when the mean vector is near the origin. In the 20 dimensional case, we numerically solve ordinary differential equation of rank 2^20 =20,148,576. We have developed software packages for the HGM. Packages based on computer algebra systems help us to solve steps (1) and (2). We have implemeted the step (3) for the Fisher-Bingham distribution, the Bingham distribution, the orthant probability, the Fisher distribution on SO(3), some of A-distributions, and the distribution function of the largest root of a Wishart matrix in the language C and/or in the system for statistics R. An implementation for the polyhedral probability is a project in progress. We find an interesting interplay with systems for polytopes in the project. References and current implementations can be found in http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hgm.html