hgm.Rhgm {hgm}

The function hgm.Rhgm performs the holonomic gradient method (HGM) for a given Pfaffian system and an initial value vector.

Description

The function hgm.Rhgm performs the holonomic gradient method (HGM) for a given Pfaffian system and an initial value vector with the deSolve package in R.

Usage

 hgm.Rhgm(th0, G0, th1, dG.fun, times=NULL, fn.params=NULL)

Arguments

th0

A d-dimensional vector which is an initial point of the parameter vector th (theta).

G0

A r-dimensional vector which is the initial value of the vector G of the normalizing constant and its derivatives.

th1

A d-dimensional vector which is the target point of th.

dG.fun

dG.fun is the “right hand sides” of the Pfaffian system. It is a d*r-dimensional array.

times

a vector; times in [0,1] at which explicit estimates for G are desired. If time = NULL, the set 0,1 is used, and only the final value is returned.

fn.params

fn.params: a list of parameters passed to the function dG.fun. If fn.params = NULL, no parameter is passed to dG.fun.

Details

The function hgm.Rhgm computes the value of a holonomic function at a given point, using HGM. This is a “Step 3” function (see the reference below), which can be used for an arbitrary input, in the HGM framework. Efficient “Step 3” functions are given for some distributions in this package.

The Pfaffian system assumed is d G_j / d th_i = (dG.fun(th, G))_i,j

The inputs of hgm.Rhgm are the initial point th0, initial value G0, final point th1, and Pfaffian system dG.fun. The output is the final value G1.

If the argument ‘times’ is specified, the function returns a matrix, where the first column denotes time, the following d-vector denotes th, and the remaining r-vector denotes G.

Value

The output is the value of G at th1. The first element of G is the normalizing constant.

Author(s)

Tomonari Sei

References

http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hgm.html

Examples

# Example 1.
# A demo program; von Mises--Fisher on S^{3-1}

G.exact = function(th){  # exact value by built-in function
  c( sinh(th[1])/th[1], cosh(th[1])/th[1] - sinh(th[1])/th[1]^2 )
}

dG.fun = function(th, G, fn.params=NULL){  # Pfaffian
  dG = array(0, c(1, 2))
  sh = G[1] * th[1]
  ch = G[2] * th[1] + G[1]
  dG[1,1] = G[2] # Pfaffian eq's
  dG[1,2] = sh/th[1] - 2*ch/th[1]^2 + 2*sh/th[1]^3
  dG
}

th0 = 0.5
th1 = 15

G0 = G.exact(th0)
G0

G1 = hgm.Rhgm(th0, G0, th1, dG.fun)  # HGM
G1

G1.exact = G.exact(th1)
G1.exact

#
# Example 2.
#
hgm.Rhgm.demo1()