[6f1cc3]: src / modules / glm / samplers / AlbertChib.cc Maximize Restore History

Download this file

AlbertChib.cc    111 lines (94 with data), 2.2 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
#include <config.h>
#include "AlbertChib.h"
#include "KS.h"
#include <graph/StochasticNode.h>
#include <graph/LinkNode.h>
#include <sampler/GraphView.h>
#include <rng/TruncatedNormal.h>
#include <rng/RNG.h>
#include <cmath>
using std::vector;
using std::string;
using std::log;
using std::exp;
using std::fabs;
//FIXME: maybe use R math library here
//Left truncated logit
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);
}
//Right truncated logit
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<GraphView const *> const &sub_views,
unsigned int chain, bool gibbs)
: BinaryGLM(view, sub_views, chain), _gibbs(gibbs)
{
}
bool AlbertChib::update(RNG *rng)
{
if (_gibbs) {
if (!updateLMGibbs(rng)) return false;
}
else {
if (!updateLM(rng)) return false;
}
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 {
return false;
//throw logic_error("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 {
return false;
//throw logic_error("Invalid child value in HolmesHeld");
}
_tau[r] = 1/sample_lambda(fabs(_z[r] - mu), rng);
break;
case BGLM_INVALID:
return false;
break;
}
}
return true;
}
string AlbertChib::name() const
{
if (_gibbs)
return "Albert-Chib-Gibbs";
else
return "Albert-Chib";
}
}