[2782e3]: src / modules / bugs / distributions / DWish.cc  Maximize  Restore  History

Download this file

220 lines (190 with data), 5.6 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
#include <config.h>
#include <rng/RNG.h>
#include <util/dim.h>
#include <util/nainf.h>
#include "lapack.h"
#include "matrix.h"
#include "DWish.h"
#include <stdexcept>
#include <cfloat>
#include <cmath>
#include <vector>
#include <JRmath.h>
using std::vector;
using std::logic_error;
using std::runtime_error;
using std::log;
#define SCALE(par) (par[0])
#define DF(par) (*par[1])
#define NROW(dims) (dims[0][0])
DWish::DWish()
: Distribution("dwish", 2, false, false)
{}
static double log_multigamma(double n, unsigned int p)
{
double y = (p * (p-1) * log(M_PI))/4;
for (unsigned int j = 0; j < p; ++j) {
y += lgamma((n-j)/2);
}
return y;
}
double DWish::logLikelihood(double const *x, unsigned int length,
vector<double const *> const &par,
vector<vector<unsigned int> > const &dims,
double const *lower, double const *upper) const
{
double const *scale = SCALE(par);
unsigned int p = NROW(dims);
double loglik = 0;
for (unsigned int i = 0; i < length; ++i) {
loglik += scale[i] * x[i];
}
loglik += DF(par) * logdet(scale, p) + (DF(par) - p - 1) * logdet(x, p);
loglik -= DF(par) * p * log(2.0) + 2 * log_multigamma(DF(par)/2, p);
return loglik/2;
}
void DWish::randomSample(double *x, int length,
double const *R, double k, int nrow,
RNG *rng)
{
/*
Generate random Wishart variable, using an algorithm proposed
by Bill Venables and originally implemented in S.
*/
if (length != nrow*nrow) {
throw logic_error("invalid length in DWish::randomSample");
}
/*
Get inverse of R. Venables' algorithm was implemented in
terms of the inverse of R, but we use a different parameterization
to preserve conjugacy.
*/
double * C = new double[length];
inverse(C, R, nrow, true);
/* Get Choleskly decomposition of C */
int info = 0;
F77_DPOTRF("U", &nrow, C, &nrow, &info);
if (info != 0) {
throw runtime_error("Failed to get Cholesky decomposition of R in dwish");
}
/* Set lower triangle of C to zero */
for (int j = 0; j < nrow; j++) {
double * C_j = &C[j*nrow]; //column j of matrix C
for (int i = j + 1; i < nrow; i++) {
C_j[i] = 0;
}
}
/* Generate square root of Wishart random variable:
- diagonal elements are square root of Chi square
- upper off-diagonal elements are normal
- lower off-diagonal elements are zero
*/
double *Z = new double[length];
for (int j = 0; j < nrow; j++) {
double *Z_j = &Z[j*nrow]; //jth column of Z
for (int i = 0; i < j; i++) {
Z_j[i] = rnorm(0, 1, rng);
}
Z_j[j] = sqrt(rchisq(k - j, rng));
for (int i = j + 1; i < nrow; i++) {
Z_j[i] = 0;
}
}
/* Transform Z with Cholesky decomposition */
double *Ztrans = new double[length];
for (int i = 0; i < nrow; i++) {
for (int j = 0; j < nrow; j++) {
double zz = 0;
for (int l = 0; l < nrow; l++) {
zz += Z[nrow * l + i] * C[nrow * j + l];
}
Ztrans[nrow * j + i] = zz;
}
}
delete [] C;
delete [] Z;
/* Now put cross-product into x */
for (int i = 0; i < nrow; i++) {
double const *Ztrans_i = &Ztrans[nrow * i];
for (int j = 0; j <= i; j++) {
double const *Ztrans_j = &Ztrans[nrow * j];
double xx = 0;
for (int l = 0; l < nrow; l++) {
xx += Ztrans_i[l] * Ztrans_j[l];
}
x[nrow * j + i] = x[nrow * i + j] = xx;
}
}
delete [] Ztrans;
}
void DWish::randomSample(double *x, unsigned int length,
vector<double const *> const &par,
vector<vector<unsigned int> > const &dims,
double const *lower, double const *upper,
RNG *rng) const
{
randomSample(x, length, SCALE(par), DF(par), NROW(dims), rng);
}
bool DWish::checkParameterDim (vector<vector<unsigned int> > const &dims) const
{
return isSquareMatrix(dims[0]) && isScalar(dims[1]);
}
vector<unsigned int>
DWish::dim(vector<vector<unsigned int> > const &dims) const
{
return dims[0];
}
bool
DWish::checkParameterValue(vector<double const *> const &par,
vector<vector<unsigned int> > const &dims) const
{
/* Check that we have sufficient degrees of freedom */
if (DF(par) < NROW(dims))
return false;
/* Check symmetry of scale matrix */
double const *scale = SCALE(par);
unsigned int nrow = NROW(dims);
for (unsigned int i = 0; i < nrow; ++i) {
for (unsigned int j = 0; j < i; ++j) {
if (fabs(scale[i + nrow*j] - scale[j + nrow*i]) > DBL_EPSILON)
return false;
}
}
/* Skipping check of positive definiteness of scale matrix */
return true;
}
void DWish::support(double *lower, double *upper, unsigned int length,
vector<double const *> const &par,
vector<vector<unsigned int> > const &dims) const
{
for (unsigned int i = 0; i < length; ++i) {
if (i % NROW(dims) == i / NROW(dims)) {
//Diagonal elements
lower[i] = 0;
}
else {
lower[i] = JAGS_NEGINF;
}
upper[i] = JAGS_POSINF;
}
}
void DWish::typicalValue(double *x, unsigned int length,
vector<double const *> const &par,
vector<vector<unsigned int> > const &dims,
double const *lower, double const *upper) const
{
/* Returns the mean as a typical value. We need to invert the
scale matrix */
inverse(x, SCALE(par), NROW(dims), true);
for (unsigned int i = 0; i < length; ++i) {
x[i] *= DF(par);
}
}
bool DWish::isSupportFixed(vector<bool> const &fixmask) const
{
return true;
}
unsigned int DWish::df(vector<vector<unsigned int> > const &dims) const
{
return dims[0][0] * (dims[0][0] + 1) / 2;
}