Digamma function’s exception, why?

megatrooperParticipant 2011-06-09 at 3:38 pm

Here is the code:

import com.numericalmethod.suanshu.analysis.uniroot.Newton;
import com.numericalmethod.suanshu.analysis.function.rn2r1.UnivariateRealFunction;
import com.numericalmethod.suanshu.analysis.function.special.Digamma;

public class Trial6{

public static void main(String[] args) {
     final double [] v = new double []{0.52392859 ,0.74728222, 0.35648189, 1.50509941, 0.09642664, 2.60895326,  0.12903240, 1.58223482, 0.08668512, 0.47263063, 0.51582472, 0.03652849, 2.33088775, 2.56285203, 0.66409021, 0.10836313, 2.11359703, 0.73836191, 0.28691570, 0.19404304, 0.22538941, 0.04399578, 1.26417153, 3.86406407, 0.12689116, 0.16198531, 0.94862916, 1.10415084, 3.39500931, 0.62369945, 3.70865055, 1.13927840, 0.12278842, 1.17936096, 2.90563559, 1.23769602};
     double [] a = new double [2];
     a=gamma_mle(v);
       Digamma digamma = new Digamma();
}

public static double [] gamma_mle (final double []y) {

UnivariateRealFunction fn = new UnivariateRealFunction() {

  @Override
           public double evaluate(double x) {
       Digamma digamma = new Digamma();
 return Math.log(x)-digamma.evaluate(x)-Math.log(sum(y)/y.length+sum(log(y))/y.length);
  }
};
        Newton newton = new Newton(fn, 1e-5);
        double alpha = newton.solve(100, 1);
        double lambda= alpha*y.length/sum(y);
        double [] vectorofpara = new double []{alpha,lambda};
        return vectorofpara;

}

public static double sum (double[] x) {
     int i,j;
     double sum=0;

    for ( j = 0; j < x.length; j++) {
        sum=sum+x[j];
    }
    return sum ;

}

public static double []  log(double[] x) {
    int j;
    for ( j = 0; j < x.length; j++) {
         x[j]= Math.log(x[j]);
    }
    return x;

}
}

The program is implemented to test the maximum likelihood estimate of gamma distribution (alpha,lambda) (supposed to be 1), given some observations (v) generated by other methods from the distribution gamma(1,1), the algorithm for finding the mle of the parameters is based on the following link from wiki : http://en.wikipedia.org/wiki/Gamma_distribution, but an error occurs due to exceptions in digamma.evaluate, how should I amend the program?

Haksun LiModerator 2011-06-09 at 6:26 pm
Hello Isaac,

I took a look at your code. There is no a problem with the
```
Digamma
```
function. However, I found a problem with your following code:
```
 
public static double [] log(double[] x) { 
 int j; 
 for ( j = 0; j < x.length; j++) { 
 x[j]= Math.log(x[j]); 
 } 
 return x; 
} 
```
I am not exactly sure what you try to do here. The problem is that it changes the input value [tt:3566cwt1]double[] x[/tt:3566cwt1].

When you call this function with [tt:3566cwt1]y[/tt:3566cwt1], you are changing the values of [tt:3566cwt1]y[/tt:3566cwt1] (the observations)!!!

It happens when you call
```
log(y)
```
in:
```
 
return Math.log(x)-digamma.evaluate(x)-Math.log(sum(y)/y.length+sum(log(y))/y.length); 
```
The first time you call
```
log(y)
```
on, e.g., [tt:3566cwt1]0.52392859[/tt:3566cwt1], it will becomes -ve. When you call
```
log(y)
```
again, it will become [tt:3566cwt1]NaN[/tt:3566cwt1].

Then, when you pass [tt:3566cwt1]NaN[/tt:3566cwt1] to the [tt:3566cwt1]Digamma[/tt:3566cwt1] function, it will throw an exception.

I suggest you modify your code to avoid changing the input values:
```
 
 private static double[] log(double[] x) { 
 int j; 
 double[] z = new double[x.length]; 
 for (j = 0; j < x.length; j++) { 
 z[j] = Math.log(x[j]); 
 } 
 return z; 
 } 
```
Alternatively, you can use the SuanShu functions: DoubleUtils.sum DoubleUtils.foreach ref: http://www.numericalmethod.com/javadoc/suanshu/com/numericalmethod/suanshu/number/DoubleUtils.html

Lastly, you should be able to debug this by walking thru you code in a debugger.

Hope this helps.