## [57223f]: src / modules / glm / samplers / AlbertChib.cc  Maximize  Restore  History

### 115 lines (98 with data), 2.5 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``` ```#include #include "AlbertChib.h" #include "KS.h" #include #include #include #include #include #include #include using std::vector; using std::string; using std::log; using std::exp; using std::fabs; //FIXME: maybe use R math library here //Random sample from a left-truncated logistic distribution static double llogit(double left, RNG *rng, double mu) { double qleft = 1/(1 + exp(mu-left)); double x = qleft + (1 - qleft) * rng->uniform(); return mu + log(x) - log(1 - x); } //Random sample from a right-truncated logistc distribution static double rlogit(double right, RNG *rng, double mu) { double qright = 1/(1 + exp(mu-right)); double x = qright * rng->uniform(); return mu + log(x) - log(1 - x); } #define CHILD(i) (_view->stochasticChildren()[i]) namespace glm { AlbertChib::AlbertChib(GraphView const *view, vector const &sub_views, unsigned int chain, bool gibbs) : BinaryGLM(view, sub_views, chain), _gibbs(gibbs), _aux_init(true) { } void AlbertChib::update(RNG *rng) { if (_aux_init) { initAuxiliary(rng); _aux_init = false; } /* Note that we must update the auxiliary variables *before* calling updateLM. This ordering is important for models with a variable design matrix (e.g. measurement error models). */ unsigned int nrow = _view->stochasticChildren().size(); double y, mu; for (unsigned int r = 0; r < nrow; ++r) { switch(_outcome[r]) { case BGLM_NORMAL: break; case BGLM_PROBIT: y = CHILD(r)->value(_chain)[0]; if (y == 1) { _z[r] = lnormal(0, rng, getMean(r)); } else if (y == 0) { _z[r] = rnormal(0, rng, getMean(r)); } else { throwLogicError("Invalid child value in HolmesHeld"); } break; case BGLM_LOGIT: y = CHILD(r)->value(_chain)[0]; mu = getMean(r); if (y == 1) { _z[r] = llogit(0, rng, mu); } else if (y == 0) { _z[r] = rlogit(0, rng, mu); } else { throwLogicError("Invalid child value in HolmesHeld"); } _tau[r] = 1/sample_lambda(fabs(_z[r] - mu), rng); break; } } if (_gibbs) { updateLMGibbs(rng); } else { updateLM(rng); } } string AlbertChib::name() const { if (_gibbs) return "Albert-Chib-Gibbs"; else return "Albert-Chib"; } } ```