Clover coverage report - baseCode - 0.2.5 Coverage timestamp: Tue Apr 12 2005 11:31:58 EDT FRAMES NO FRAMES |
file stats: | LOC: | 159 |  | Methods: | 2
| NCLOC: | 43 |  | Classes: | 1
30 day Evaluation Version distributed via the Maven Jar Repository. Clover is not free. You have 30 days to evaluate it. Please visit http://www.thecortex.net/clover to obtain a licensed version of Clover
Source file | Conditionals | Statements | Methods | TOTAL
Wishart.java | 70% | 86.4% | 100% | 82.4% |

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 | 1 |
        public Wishart( double df, DoubleMatrix2D covariance, RandomEngine randomGenerator ) {
44 | 1 |
           this.s = covariance.columns();
45 | 0 |
           if ( s != covariance.rows() ) throw new IllegalArgumentException( "Covariance matrix must be square" );
46 | 0 |
           if ( df <= s - 1 ) throw new IllegalArgumentException( "df must be greater than s - 1" );
47 | 0 |
           if ( randomGenerator == null ) throw new IllegalArgumentException( "Null random number generator" );
48 |  |
49 | 1 |
           this.r = randomGenerator;
50 | 1 |
           this.df = df;
51 | 1 |
           this.cov = covariance;
52 |  |
53 | 1 |
           rgamma = new Gamma( 1, 1, r );
54 | 1 |
           rnorm = new Normal( s * ( s - 1.0 ) / 2.0, 1.0, r );
55 |  |
56 | 1 |
           CholeskyDecomposition c = new CholeskyDecomposition( covariance );
57 | 1 |
           chol = a.transpose( c.getL() ); // returns lower triangle so we transpose to make upper triangular.
58 | 1 |
           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 | 1 |
        public DoubleMatrix2D nextDoubleMatrix() {
118 | 1 |
           mat.assign( 0.0 );
119 |  |
120 |  |
           // fill in diagonal with random gamma deviates, upper triangle with random normal deviates.
121 | 1 |
           for ( int i = 0; i < s; i++ ) {
122 | 2 |
              mat.setQuick( i, i, Math.sqrt( 2 * rgamma.nextDouble( s, ( df + 1.0 - i ) / 2.0 ) ) );
123 | 2 |
              for ( int j = i + 1; j < s; j++ ) {
124 | 1 |
                 mat.setQuick( i, j, rnorm.nextDouble() );
125 |  |
              }
126 |  |
           }
127 |  |
128 | 1 |
           mat = a.mult( mat, chol );
129 | 1 |
           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 |  |
     }