Evaluating the normalizing constant of the Fisher-Bingham distribution by HGM and solving MLE problems by HGD

Programs and data used in the paper
T. Koyama, H. Nakayama, K. Nishiyama, N. Takayama, Holonomic Gradient Descent for the Fisher-Bingham Distribution on the n-dimensional Sphere, arxiv:1201.3239
are here (fb2.zip) in the zip format.
In order to compile and execute them, you need the Gnu Scientific Library (GSL) and R software .
( Risa/Asir , contained in the openxm package, is optimal to execute .rr programs.)
Sample data used in the paper are in the folder fb2/ex in the zip file.

An example on S^3 (the problem s3_e1 in the paper)

  <-- unzip fb2.zip
cd fb2
make all
R        <-- start R
> source("s_test.r");
> call_s3_e1();
Read 17 items            <-- 1st step by ko-initial and optim()
[1] 26.21370 -3.76030  1.70057  3.39803  5.16930  5.25944  4.58366  8.34724
[1] 23.5735
Read 17 items
[1] 26.72370 -3.91633  1.72777  3.44281  5.24418  5.48767  4.65147  8.45515
[1] 23.68049
Read 17 items
[1] 26.18900 -3.73436  1.65283  3.39216  5.16862  5.25063  4.57549  8.34269
[1] 23.58164

.... 

Read 17 items
[1] 46.66620  1.89242  3.51563  9.26612 -9.51110  6.50483  8.77901 15.13790
[1] 10.08021

         <-- The output of optim() is in the file s3_e1_pos_r.

FR method               <-- 2nd step by HGD (C program s3_e1)
pos[0] = -2.108446
pos[1] = -0.814766
pos[2] = -3.796403
pos[3] = 0.753867
.... snip
1, -2.10603, -0.817113, -3.79343, 0.748836, -1.15559, -0.125705, 4.16824, 1.01524, -0.543168, 1.05994, 1.18925, 1.28324, 0.893618, -1.14147, 10.0525
points = [-2.10603, -0.817113, -3.79343, 0.748836, -1.15559, -0.125705, 4.16824, 1.01524, -0.543168, 1.05994, 1.18925, 1.28324, 0.893618, -1.14147]
values = [10.0525, 0.47601, 0.703612, 1.82723, -2.11336, 1.41809, 1.87452, 3.20026]
grad : -0.142887 0.141284 -0.178817 0.299841 0.013521 0.012832 -0.129449 -0.020050 -0.000856 0.149416 -0.169592 -0.107367 0.291957 -0.108757 
norm(grad): 0.582620
2, -2.10121, -0.821808, -3.78748, 0.738776, -1.15595, -0.126147, 4.17259, 1.01583, -0.542996, 1.05488, 1.19527, 1.28701, 0.883296, -1.13746, 10.0413
points = [-2.10121, -0.821808, -3.78748, 0.738776, -1.15595, -0.126147, 4.17259, 1.01583, -0.542996, 1.05488, 1.19527, 1.28701, 0.883296, -1.13746]
values = [10.0413, 0.497838, 0.714846, 1.7884, -2.08791, 1.42095, 1.87727, 3.19244]
grad : -0.138283 0.141407 -0.178792 0.294271 0.018362 0.011912 -0.126660 -0.024261 0.007705 0.144182 -0.147041 -0.095226 0.254850 -0.085553 
norm(grad): 0.546439

.... snip

norm(grad): 0.001030
Minimum found        <--
73, -1.22533, -1.10699, -2.63975, -0.0812375, -1.85274, 0.00702786, 5.35051, 1.31406, -0.929157, 0.577568, 0.887881, 1.65851, 0.593832, -1.31242, 9.58807
points = [-1.22533, -1.10699, -2.63975, -0.0812375, -1.85274, 0.00702786, 5.35051, 1.31406, -0.929157, 0.577568, 0.887881, 1.65851, 0.593832, -1.31242] <-- the optimum point
values = [9.58807, 0.615897, 0.773816, 1.46424, -1.91159, 1.48891, 1.77486, 3.07148] <-- the optimum value
grad : 0.000061 0.000164 -0.000315 -0.000037 -0.000148 -0.000176 -0.000431 -0.000033 -0.000299 0.000121 0.000125 0.000307 -0.000101 0.000390 
norm(grad): 0.000864
search_min :
[9.58807, 0.615897, 0.773816, 1.46424, -1.91159, 1.48891, 1.77486, 3.07148]
time(search_min) : 3.240000
... snip 
>