Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

[e23e50]: src / modules / glm / samplers / HolmesHeld.cc Maximize Restore History

Download this file

HolmesHeld.cc    219 lines (177 with data), 6.0 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 "HolmesHeld.h"
#include "KS.h"
#include "Outcome.h"
#include <graph/StochasticNode.h>
#include <graph/LinkNode.h>
#include <sampler/GraphView.h>
#include <rng/TruncatedNormal.h>
#include <rng/RNG.h>
#include <module/ModuleError.h>
extern "C" {
#include <cs.h>
}
#include <cmath>
// Regularization penalty for precision
#define REG_PENALTY 0.001
using std::vector;
using std::string;
using std::sqrt;
extern cholmod_common *glm_wk;
static cs shallow_copy(cholmod_sparse *s)
{
//Copy a cholmod_sparse struct to a cs struct without allocating
//any memory
cs c;
c.nzmax = s->nzmax;
c.n = s->nrow;
c.m = s->ncol;
c.p = static_cast<int*>(s->p);
c.i = static_cast<int*>(s->i);
c.x = static_cast<double*>(s->x);
c.nz = -1;
return c;
}
namespace jags {
namespace glm {
HolmesHeld::HolmesHeld(GraphView const *view,
vector<GraphView const *> const &sub_views,
vector<Outcome *> const &outcomes,
unsigned int chain)
: GLMMethod(view, sub_views, outcomes, chain), _aux_init(true)
{
}
void HolmesHeld::updateAuxiliary(cholmod_dense *w,
cholmod_factor *N, RNG *rng)
{
/*
In the parent GLMMethod class, the posterior precision of
the regression parameters is represented by the matrix "A";
the posterior mean "mu" solves A %*% mu = b.
In this call, "N" holds the factorization of P %*% A %*% t(P),
where P is a permutation matrix. The factorization takes
the form L %*% D %*% t(L), where D is diagonal and L is
a lower triangular matrix. The parameter "w" solves
L %*% w = P %*% b
IMPORTANT NOTE: mu, b use a parameterization in which the
current value of the regression parameters is taken as the
origin. This requires us to adjust the calculations using
the current value of the linear predictor (mu_r).
*/
vector<StochasticNode *> const &schildren =
_view->stochasticChildren();
unsigned int nrow = schildren.size();
//Transpose and permute the design matrix
cholmod_sparse *t_x = cholmod_transpose(_x, 1, glm_wk);
int *fperm = static_cast<int *>(_factor->Perm);
cholmod_sparse *Pt_x = cholmod_submatrix(t_x, fperm, t_x->nrow,
0, -1, 1, 1, glm_wk);
cholmod_free_sparse(&t_x, glm_wk);
//Turn the factor into a matrix. Since this modifies the factor
//we need to take a copy first
cholmod_factor *f = cholmod_copy_factor(_factor, glm_wk);
cholmod_sparse *L = cholmod_factor_to_sparse(f, glm_wk);
if (!L->packed || !L->sorted) {
throwLogicError("Cholesky factor is not packed or not sorted");
}
cholmod_free_factor(&f, glm_wk);
unsigned int ncol = L->ncol;
vector<double> d(ncol, 1);
if (!_factor->is_ll) {
//Extract D from the diagonal of L and set the diagonal to 1
int *Lp = static_cast<int*>(L->p);
double *Lx = static_cast<double*>(L->x);
for (unsigned int r = 0; r < ncol; ++r) {
d[r] = Lx[Lp[r]];
Lx[Lp[r]] = 1;
}
}
// Now make shallow copies of L and Pt_x so we can pass them
// to cs_sparse from the CSparse library
cs cs_L = shallow_copy(L);
cs cs_Ptx = shallow_copy(Pt_x);
double *ur = new double[ncol];
int *xi = new int[2*ncol]; //Stack
double *wx = static_cast<double *>(w->x);
for (unsigned int r = 0; r < nrow; ++r) {
// In a heterogeneous GLM, we may have some normal
// outcomes as well as binary outcomes. These can be
// skipped as there is no need for auxiliary variables
/* FIXME
No way to skip in new API
if (_outcomes[r] == BGLM_NORMAL)
continue;
*/
int top = cs_spsolve(&cs_L, &cs_Ptx, r, xi, ur, 0, 1);
double mu_r = _outcomes[r]->mean(); // See IMPORTANT NOTE above
double tau_r = _outcomes[r]->precision();
//Calculate mean and precision of z[r] conditional
//on z[s] for s != r
double zr_mean = 0;
double Hr = 0; // Diagonal of hat matrix
if (_factor->is_ll) {
for (unsigned int j = top; j < ncol; ++j) {
zr_mean += ur[xi[j]] * wx[xi[j]];
Hr += ur[xi[j]] * ur[xi[j]];
}
}
else {
for (unsigned int j = top; j < ncol; ++j) {
zr_mean += ur[xi[j]] * wx[xi[j]] / d[xi[j]];
Hr += ur[xi[j]] * ur[xi[j]] / d[xi[j]];
}
}
Hr *= tau_r;
if (1 - Hr <= 0) {
// Theoretically this should never happen, but numerical instability may cause it
StochasticNode const *snode = _view->stochasticChildren()[r];
throwNodeError(snode, "Highly influential outcome variable in Holmes-Held update method.");
}
zr_mean -= Hr * (_outcomes[r]->value()_z[r] - mu_r);
zr_mean /= (1 - Hr);
double zr_prec = (1 - Hr) * tau_r;
double yr = schildren[r]->value(_chain)[0];
double zold = _z[r];
if (yr == 1) {
_z[r] = lnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
}
else if (yr == 0) {
_z[r] = rnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
}
else {
throwLogicError("Invalid child value in HolmesHeld");
}
//Add new contribution of row r back to b
double zdelta = (_z[r] - zold) * tau_r;
for (unsigned int j = top; j < ncol; ++j) {
wx[xi[j]] += ur[xi[j]] * zdelta;
}
}
//Free workspace
delete [] ur;
delete [] xi;
cholmod_free_sparse(&Pt_x, glm_wk);
cholmod_free_sparse(&L, glm_wk);
}
void HolmesHeld::update(RNG *rng)
{
if (_aux_init) {
for (unsigned int r = 0; r < _outcomes.size(); ++r) {
_outcomes[r]->update(rng);
}
_aux_init = false;
}
/* Well, this is awkward. Nothing in the new API to update mixture parameters separate
from the other auxiliary variables
FIXME
for (unsigned int r = 0; r < _tau.size(); ++r)
{
if (_outcome[r] == BGLM_LOGIT) {
double delta = fabs(_outcomes[r]->value() - _outcomes[r]->mean());
_tau[r] = REG_PENALTY + 1/sample_lambda(delta, rng);
}
}
*/
updateLM(rng);
}
}}