View Javadoc

    1   package baseCode.math.distribution;
    2
    3   import cern.colt.list.DoubleArrayList;
    4   import cern.colt.matrix.DoubleMatrix2D;
    5   import cern.colt.matrix.impl.DenseDoubleMatrix2D;
    6   import cern.colt.matrix.linalg.Algebra;
    7   import cern.jet.random.Beta;
    8   import cern.jet.random.engine.RandomEngine;
    9   import cern.jet.stat.Descriptive;
    10
    11  /***
    12   * <hr>
    13   * <p>
    14   * Copyright (c) 2004 Columbia University
    15   *
    16   * @author pavlidis
    17   * @version $Id: Dirichlet.java,v 1.1 2005/03/21 18:01:03 pavlidis Exp $
    18   */
    19  public class Dirichlet {
    20
    21     private Algebra a = new Algebra();
    22     private RandomEngine r;
    23
    24     private DoubleArrayList p;
    25     private Beta rbeta;
    26
    27     public Dirichlet( DoubleArrayList p, RandomEngine randomGenerator ) {
    28        if ( randomGenerator == null ) throw new IllegalArgumentException( "Null random number generator" );
    29        this.r = randomGenerator;
    30        this.p = p;
    31        rbeta = new Beta( 1, 1, r );
    32     }
    33
    34     public double nextDouble() {
    35
    36        return 0;
    37     }
    38
    39     public DoubleArrayList draws( int n ) {
    40
    41        DoubleMatrix2D mat = new DenseDoubleMatrix2D( n, p.size() );
    42        double psum = Descriptive.sum( (DoubleArrayList)p.partFromTo( 1, p.size() - 1 ) );
    43
    44        for ( int i = 0; i < mat.rows(); i++ ) {
    45           mat.setQuick( 0, i, rbeta.nextDouble( p.getQuick( 0 ), psum ) );
    46        }
    47
    48        for (int i = 1; i < p.size(); i++) {
    49           for ( int j = 0; j < mat.rows(); j++ ) {
    50              mat.setQuick( i, j , rbeta.nextDouble( p.getQuick( 0 ), psum ) );
    51           }
    52        }
    53
    54        return null;
    55     }
    56  }
    57  //rdirichlet <- function ( n, p ) {
    58  //   # return n random samples from a Dirichlet distribution with parameter p
    59  //     if ( !is.vector(n, "numeric")
    60  //          | length(n) != 1
    61  //          | !is.vector(p, "numeric")
    62  //        ) { stop("error in call to rdirichlet") }
    63  //     mat <- matrix ( NA, n, length(p) )
    64  //     mat[,1] <- rbeta ( n, p[1], sum(p[-1]) )
    65  //
    66  //   for ( i in 2:(length(p)-1) ) {
    67  //
    68  //     mat[,i] <- ( rbeta ( n, p[i], sum(p[(i+1):length(p)]) )
    69  //                     * ( 1 - apply ((mat[,1:(i-1),drop=F]), 1, sum) )
    70  //                  )
    71  //     }
    72  //     mat[,length(p)] <- 1 - apply ( (mat[,-length(p),drop=F]), 1, sum )
    73  //
    74  //     return ( mat )
    75  //   }