[4c0f6f]: src / modules / bugs / samplers / ConjugateWishart.cc Maximize Restore History

Download this file

ConjugateWishart.cc    154 lines (126 with data), 3.8 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
#include <config.h>
#include "ConjugateWishart.h"
#include "DWish.h"
#include <rng/RNG.h>
#include <graph/LogicalNode.h>
#include <graph/StochasticNode.h>
#include <graph/MixtureNode.h>
#include <sampler/Linear.h>
#include <sampler/GraphView.h>
#include <set>
#include <vector>
#include <cmath>
#include <algorithm>
#include <JRmath.h>
using std::vector;
using std::set;
using std::sqrt;
using std::string;
using std::copy;
static inline double getPrecision0(StochasticNode const *snode,
unsigned int chain)
{
//Returns the first element of the precision matrix for a node
//with a multivariate normal distribution.
return snode->parents()[1]->value(chain)[0];
}
namespace bugs {
bool ConjugateWishart::canSample(StochasticNode *snode, Graph const &graph)
{
if (getDist(snode) != WISH)
return false;
if (isBounded(snode))
return false;
GraphView gv(snode, graph);
vector<StochasticNode const*> const &schild = gv.stochasticChildren();
// Check stochastic children
for (unsigned int i = 0; i < schild.size(); ++i) {
if (isBounded(schild[i])) {
return false; //Bounded
}
switch(getDist(schild[i])) {
case MNORM:
if (gv.isDependent(schild[i]->parents()[0])) {
return false; //mean parameter depends on snode
}
break;
default:
return false;
}
}
vector<DeterministicNode*> const &dchild = gv.deterministicChildren();
if (!dchild.empty()) {
// Deterministic children must be scale functions
if (!checkScale(&gv, false)) {
return false;
}
// Only mixture nodes are allowed. If we allowed arbitrary
// functions, the complexity would be O(nrow^4).
for (unsigned int i = 0; i < dchild.size(); ++i) {
if (!isMixture(dchild[i]))
return false;
}
}
return true;
}
ConjugateWishart::ConjugateWishart(GraphView const *gv)
: ConjugateMethod(gv)
{}
void
ConjugateWishart::update(unsigned int chain, RNG *rng) const
{
vector<StochasticNode const*> const &stoch_children =
_gv->stochasticChildren();
unsigned int nchildren = stoch_children.size();
vector<Node const *> const &param = _gv->nodes()[0]->parents();
double df = *param[1]->value(chain);
double const *Rprior = param[0]->value(chain);
int nrow = param[0]->dim()[0];
int N = nrow * nrow;
vector<double> R(N);
copy(Rprior, Rprior + N, R.begin());
//Logical mask to determine which stochastic children are active.
vector<bool> active(nchildren, true);
if (!_gv->deterministicChildren().empty()) {
//Mixure model
//Save first element of precision matrix for each child
vector<double> precision0(nchildren);
for (unsigned int i = 0; i < nchildren; ++i) {
precision0[i] = getPrecision0(stoch_children[i], chain);
}
//Double the current value
double const *x = _gv->nodes()[0]->value(chain);
vector<double> x2(N);
for (int j = 0; j < N; ++j) {
x2[j] = 2 * x[j];
}
_gv->setValue(x2, chain);
//See if precision matrix has changed
for (unsigned int i = 0; i < nchildren; ++i) {
if (getPrecision0(stoch_children[i], chain) == precision0[i]) {
active[i] = false; //not active
}
}
}
for (unsigned int i = 0; i < nchildren; ++i) {
if (active[i]) {
StochasticNode const *schild = stoch_children[i];
double const *Y = schild->value(chain);
double const *mu = schild->parents()[0]->value(chain);
for (int j = 0; j < nrow; j++) {
for (int k = 0; k < nrow; k++) {
R[j*nrow + k] += (Y[j] - mu[j]) * (Y[k] - mu[k]);
}
}
df += 1;
}
}
vector<double> xnew(N);
DWish::randomSample(&xnew[0], N, &R[0], df, nrow, rng);
_gv->setValue(xnew, chain);
}
string ConjugateWishart::name() const
{
return "ConjugateWishart";
}
}