Diff of /src/modules/glm/samplers/HolmesHeld.cc [01a998] .. [081fbe] Maximize Restore

  Switch to unified view

a/src/modules/glm/samplers/HolmesHeld.cc b/src/modules/glm/samplers/HolmesHeld.cc
...
...
13
extern "C" {
13
extern "C" {
14
#include <cs.h>
14
#include <cs.h>
15
}
15
}
16
16
17
#include <cmath>
17
#include <cmath>
18
19
// Regularization penalty for precision 
20
#define REG_PENALTY 0.001
18
21
19
using std::vector;
22
using std::vector;
20
using std::string;
23
using std::string;
21
using std::sqrt;
24
using std::sqrt;
22
25
...
...
49
52
50
    void HolmesHeld::updateAuxiliary(cholmod_dense *w, 
53
    void HolmesHeld::updateAuxiliary(cholmod_dense *w, 
51
                     cholmod_factor *N, RNG *rng)
54
                     cholmod_factor *N, RNG *rng)
52
    {
55
    {
53
    /* 
56
    /* 
54
       In the parent GLMMethod class, the posterior precision is
57
       In the parent GLMMethod class, the posterior precision of
55
     represented by the matrix "A"; the posterior mean "mu"
58
           the regression parameters is represented by the matrix "A";
56
     solves A %*% mu = b.
59
           the posterior mean "mu" solves A %*% mu = b.
57
       
60
       
58
     In this call, "w" solves A %*% w = P %*% b and "N" holds
61
           In this call, "N" holds the factorization of P %*% A %*% t(P), 
59
     the Cholesky decomposition of P %*% A %*% t(P). 
62
           where P is a permutation matrix.  The factorization takes
63
           the form L %*% D %*% t(L), where D is diagonal and L is
64
           a lower triangular matrix. The parameter "w" solves
65
           L %*% w = P %*% b
66
     
67
     IMPORTANT NOTE: mu, b use a parameterization in which the
68
     current value of the regressions parameters is taken as the
69
     origin.  This requires us to adjust the calculations using
70
     the current value of the linear predictor (mu_r).
60
    */
71
    */
61
72
62
    vector<StochasticNode const *> const &schildren = 
73
    vector<StochasticNode const *> const &schildren = 
63
        _view->stochasticChildren();
74
        _view->stochasticChildren();
64
75
...
...
102
        int *xi = new int[2*ncol]; //Stack
113
        int *xi = new int[2*ncol]; //Stack
103
    
114
    
104
    double *wx = static_cast<double *>(w->x);
115
    double *wx = static_cast<double *>(w->x);
105
    
116
    
106
    for (unsigned int r = 0; r < nrow; ++r) {
117
    for (unsigned int r = 0; r < nrow; ++r) {
107
      
118
119
      // In a heterogeneous GLM, we may have some normal
120
      // outcomes as well as binary outcomes. These can be
121
      // skipped as there is no need for auxiliary variables
108
        if (_outcome[r] != BGLM_NORMAL) {
122
        if (_outcome[r] == BGLM_NORMAL)
123
      continue;
109
        
124
        
110
      int top = cs_spsolve(&cs_L, &cs_Ptx, r, xi, ur, 0, 1);
125
        int top = cs_spsolve(&cs_L, &cs_Ptx, r, xi, ur, 0, 1);
111
        
126
        
112
      double mu_r = getMean(r);
127
      double mu_r = getMean(r); // See IMPORTANT NOTE above
113
      double tau_r = getPrecision(r);
128
        double tau_r = getPrecision(r);
114
129
      
115
      //Calculate mean and precision of z[r] conditional
130
        //Calculate mean and precision of z[r] conditional
116
      //on z[s] for s != r
131
        //on z[s] for s != r
117
      double zr_mean = 0;
132
        double zr_mean = 0;
118
      double Hr = 0; // 
133
      double Hr = 0; // Diagonal of hat matrix
119
134
      
120
      if (_factor->is_ll) {
135
        if (_factor->is_ll) {
121
            for (unsigned int j = top; j < ncol; ++j) {
136
        for (unsigned int j = top; j < ncol; ++j) {
122
          zr_mean  += ur[xi[j]] * wx[xi[j]]; 
137
            zr_mean  += ur[xi[j]] * wx[xi[j]]; 
123
          Hr  += ur[xi[j]] * ur[xi[j]];
138
            Hr  += ur[xi[j]] * ur[xi[j]];
124
          }
125
        }
139
        }
140
      }
126
      else {
141
        else {
127
            for (unsigned int j = top; j < ncol; ++j) {
142
        for (unsigned int j = top; j < ncol; ++j) {
128
          zr_mean  += ur[xi[j]] * wx[xi[j]] / d[xi[j]]; 
143
            zr_mean  += ur[xi[j]] * wx[xi[j]] / d[xi[j]]; 
129
          Hr  += ur[xi[j]] * ur[xi[j]] / d[xi[j]];
144
            Hr  += ur[xi[j]] * ur[xi[j]] / d[xi[j]];
130
          }
131
        }
145
        }
146
      }
132
      Hr *= tau_r;
147
        Hr *= tau_r;
148
      if (1 - Hr <= 0) {
149
      // Theoretically this should never happen, but numerical instability may cause it
150
      StochasticNode const *snode = _view->stochasticChildren()[r];
151
      throwNodeError(snode, "Highly influential outcome variable in Holmes-Held update method.");
152
      }
133
      zr_mean -= Hr * (_z[r] - mu_r);
153
        zr_mean -= Hr * (_z[r] - mu_r);
134
      zr_mean /= (1 - Hr);
154
        zr_mean /= (1 - Hr);
135
      double zr_prec = (1 - Hr) * tau_r;
155
        double zr_prec = (1 - Hr) * tau_r;
136
      
156
      
137
      if (zr_prec <= 0) {
138
          throwRuntimeError("Invalid precision in Holmes-Held update method.\nThis is a known bug and we are working on it.\nPlease bear with us");
139
      }
140
141
      double yr = schildren[r]->value(_chain)[0];
157
        double yr = schildren[r]->value(_chain)[0];
142
      double zold = _z[r];
158
        double zold = _z[r];
143
      if (yr == 1) {
159
        if (yr == 1) {
144
            _z[r] = lnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
160
        _z[r] = lnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
145
      }
161
      }
146
      else if (yr == 0) {
162
        else if (yr == 0) {
147
            _z[r] = rnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
163
        _z[r] = rnormal(0, rng, mu_r + zr_mean, 1/sqrt(zr_prec));
148
      }
164
      }
149
      else {
165
        else {
150
            throwLogicError("Invalid child value in HolmesHeld");
166
        throwLogicError("Invalid child value in HolmesHeld");
151
      }
167
      }
152
168
      
153
      //Add new contribution of row r back to b
169
        //Add new contribution of row r back to b
154
      double zdelta = (_z[r] - zold) * tau_r;
170
        double zdelta = (_z[r] - zold) * tau_r;
155
      for (unsigned int j = top; j < ncol; ++j) {
171
        for (unsigned int j = top; j < ncol; ++j) {
156
            wx[xi[j]] += ur[xi[j]] * zdelta; 
172
        wx[xi[j]] += ur[xi[j]] * zdelta; 
157
      }
158
        }
173
        }
159
    }
174
    }
160
        
175
        
161
    //Free workspace
176
    //Free workspace
162
    delete [] ur;
177
    delete [] ur;
163
    delete [] xi;
178
    delete [] xi;
164
    
179
    
165
  //cholmod_free_sparse(&u, glm_wk);
166
    cholmod_free_sparse(&Pt_x, glm_wk);
180
    cholmod_free_sparse(&Pt_x, glm_wk);
167
    cholmod_free_sparse(&L, glm_wk);
181
    cholmod_free_sparse(&L, glm_wk);
168
    }
182
    }
169
    
183
    
170
    void HolmesHeld::update(RNG *rng)
184
    void HolmesHeld::update(RNG *rng)
...
...
172
    if (_aux_init) {
186
    if (_aux_init) {
173
        initAuxiliary(rng);
187
        initAuxiliary(rng);
174
        _aux_init = false;
188
        _aux_init = false;
175
    }
189
    }
176
190
177
  /*
178
    Update the auxiliary variables *before* calling
179
    updateLM. This ordering is important for models with a
180
    variable design matrix (e.g.  measurement error models).
181
  */
182
    for (unsigned int r = 0; r < _tau.size(); ++r)
191
    for (unsigned int r = 0; r < _tau.size(); ++r)
183
    {
192
    {
184
        if (_outcome[r] == BGLM_LOGIT) {
193
        if (_outcome[r] == BGLM_LOGIT) {
185
        double delta = fabs(getValue(r) - getMean(r));
194
        double delta = fabs(getValue(r) - getMean(r));
186
        _tau[r] = 1/sample_lambda(delta, rng);
195
        _tau[r] = REG_PENALTY + 1/sample_lambda(delta, rng);
187
        }
196
        }
188
    }
197
    }
189
198
190
    updateLM(rng);
199
    updateLM(rng);
191
    }
200
    }