Login
Login

Digamma function’s exception, why?

Home 21090308 Forums Numerical Method Analysis Digamma function’s exception, why?

  • This topic is empty.
Viewing 1 post (of 1 total)
  • Author
    Posts
  • #1871
    megatrooper
    Member

    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?

Viewing 1 post (of 1 total)
  • You must be logged in to reply to this topic.