View Javadoc

    1   package baseCode.math.distribution;
    2
    3   import cern.colt.matrix.DoubleMatrix2D;
    4   import cern.colt.matrix.impl.DenseDoubleMatrix2D;
    5   import cern.colt.matrix.linalg.Algebra;
    6   import cern.colt.matrix.linalg.CholeskyDecomposition;
    7   import cern.jet.random.Gamma;
    8   import cern.jet.random.Normal;
    9   import cern.jet.random.engine.RandomEngine;
    10
    11  /***
    12   * Wishart distribution, used to simulate covariance matrices.
    13   * <p>
    14   * Based on method in Odell and Feiveson JASA 1966 p.199-203
    15   * <p>
    16   * The interface is modeled after ContinuousDistribution from colt, which unfortunately is designed only for univariate
    17   * distributions.
    18   * <hr>
    19   * <p>
    20   * Copyright (c) 2004 Columbia University
    21   *
    22   * @author pavlidis
    23   * @version $Id: Wishart.java,v 1.1 2005/01/16 03:51:11 pavlidis Exp $
    24   */
    25  public class Wishart {
    26
    27     private Gamma rgamma;
    28     int s; // dimension of matrix
    29     double df; // degrees of freedom
    30     DoubleMatrix2D cov; // input covariance matrix
    31     DoubleMatrix2D chol; // cholesky decomposition of the covariance matrix.
    32     private DoubleMatrix2D mat;
    33     private Normal rnorm;
    34     private Algebra a = new Algebra();
    35     private RandomEngine r;
    36
    37     /***
    38      * @param s
    39      * @param df
    40      * @param covariance
    41      * @param randomGenerator
    42      */
    43     public Wishart( double df, DoubleMatrix2D covariance, RandomEngine randomGenerator ) {
    44        this.s = covariance.columns();
    45        if ( s != covariance.rows() ) throw new IllegalArgumentException( "Covariance matrix must be square" );
    46        if ( df <= s - 1 ) throw new IllegalArgumentException( "df must be greater than s - 1" );
    47        if ( randomGenerator == null ) throw new IllegalArgumentException( "Null random number generator" );
    48
    49        this.r = randomGenerator;
    50        this.df = df;
    51        this.cov = covariance;
    52
    53        rgamma = new Gamma( 1, 1, r );
    54        rnorm = new Normal( s * ( s - 1.0 ) / 2.0, 1.0, r );
    55
    56        CholeskyDecomposition c = new CholeskyDecomposition( covariance );
    57        chol = a.transpose( c.getL() ); // returns lower triangle so we transpose to make upper triangular.
    58        mat = new DenseDoubleMatrix2D( this.s, this.s );
    59     }
    60
    61     /***
    62      * Based on R code from Francesca Dominici, <a
    63      * href="http://www.biostat.jhsph.edu/~fdominic/teaching/BM/bm.html">http://www.biostat.jhsph.edu/~fdominic/teaching/BM/bm.html
    64      * </a>
    65      * <p>
    66      * Returns
    67      *
    68      * <pre>
    69      *
    70      *
    71      *
    72      *
    73      *
    74      *
    75      *                 w=(RU)'RU
    76      *
    77      *
    78      *
    79      *
    80      *
    81      *
    82      * </pre>
    83      *
    84      * where
    85      *
    86      * <pre>
    87      *
    88      *
    89      *
    90      *
    91      *
    92      *
    93      *
    94      *
    95      *         Cov=U'U (U is upper triang)
    96      *
    97      *
    98      *
    99      * </pre>
    100     *
    101     * and where upper-tri R is
    102     *
    103     * <pre>
    104     *
    105     *
    106     *                                         R_ij&tilde;N(0,1), i&lt;j ; (R_ii)&circ;2&tilde;Chisq(nu-s+i)
    107     *
    108     *
    109     *
    110     * </pre>
    111     *
    112     * @param s
    113     * @param nu
    114     * @param covariance
    115     * @return
    116     */
    117    public DoubleMatrix2D nextDoubleMatrix() {
    118       mat.assign( 0.0 );
    119
    120       // fill in diagonal with random gamma deviates, upper triangle with random normal deviates.
    121       for ( int i = 0; i < s; i++ ) {
    122          mat.setQuick( i, i, Math.sqrt( 2 * rgamma.nextDouble( s, ( df + 1.0 - i ) / 2.0 ) ) );
    123          for ( int j = i + 1; j < s; j++ ) {
    124             mat.setQuick( i, j, rnorm.nextDouble() );
    125          }
    126       }
    127
    128       mat = a.mult( mat, chol );
    129       return a.mult( a.transpose( mat ).copy(), mat );
    130
    131    }
    132
    133    //   #GENERATE WISHART------------------------------------------------------------
    134    //   "rwish" <- function(s,nu,Cov)
    135    //   {
    136    //   #sxs Wishart matrix, nu degree of freedom, var/covar Cov based on
    137    //   #P.L.Odell & A.H. Feiveson(JASA 1966 p.199-203). Returns w=(RU)'RU
    138    //   #where Cov=U'U (U is upper triang) and where upper-tri R is
    139    //   # R_ij~N(0,1), i<j ; (R_ii)^2~Chisq(nu-s+i)
    140    //   if (nu<=s-1) stop ("Wishart algorithm requires nu>s-1")
    141    //   R<- diag(sqrt(2*rgamma(s,(nu+1 - 1:s)/2)))
    142    //   R[outer(1:s, 1:s, "<")] <- rnorm (s*(s-1)/2)
    143    //   R <- R%*% chol(Cov)
    144    //   return(t(R)%*%R)
    145    //   }
    146    //   #GENERATE INVERSE WISHART----------------------------------------------------
    147    //   "riwish" <- function(s,df,Prec)
    148    //   {
    149    //   #sxs Inverse Wishart matrix, df degree of freedom, precision matrix
    150    //   #Prec. Distribution of W^{-1} for Wishart W with nu=df+s-1 degree of
    151    //   # freedoom, covar martix Prec^{-1}.
    152    //   # NOTE mean of riwish is proportional to Prec
    153    //   if (df<=0) stop ("Inverse Wishart algorithm requires df>0")
    154    //   R <- diag(sqrt(2*rgamma(s,(df + s - 1:s)/2)))
    155    //   R[outer(1:s, 1:s, "<")] <- rnorm (s*(s-1)/2)
    156    //   S <- t(solve(R))%*% chol(Prec)
    157    //   return(t(S)%*%S)
    158
    159 }