#include #include #include #include #include #include #include #include #include using std::logic_error; using std::string; FiniteSampler::FiniteSampler(StochasticNode *node, Graph const &graph) : GibbsSampler(node, graph) { if (!canSample(node, graph)) throw logic_error("Attempt to construct invalid FiniteSampler"); } FiniteSampler::~FiniteSampler() {} void FiniteSampler::update() { Distribution const *dist = node()->distribution(); long lower = static_cast(dist->lowerSupport(0, node()->parameters())); long upper = static_cast(dist->upperSupport(0, node()->parameters())); long size = upper - lower + 1; double *lik = new double[size]; double liksum = 0.0; for (long i = 0; i < size; i++) { double ivalue = lower + i; setValue(&ivalue, 1); lik[i] = exp(logFullConditional()); liksum += lik[i]; } /* Sample */ double urand = runif(0.0, liksum); long i; liksum = 0.0; for (i = 0; i < size - 1; i++) { liksum += lik[i]; if (liksum > urand) { break; } } double ivalue = lower + i; setValue(&ivalue, 1); delete [] lik; } bool FiniteSampler::canSample(StochasticNode const *node, Graph const &graph) { //Node must be scalar with discrete-valued distribution of full rank Distribution const *dist = node->distribution(); if (!dist->isDiscreteValued()) return false; if (node->data.length() != 1) return false; if (dist->df(node->parameters()) != 1) return false; //Distribution cannot be unbounded if (dist->upperSupport(0, node->parameters()) == DBL_MAX || dist->lowerSupport(0, node->parameters()) == -DBL_MAX) { return false; } else { double n = dist->upperSupport(0, node->parameters()) - dist->lowerSupport(0, node->parameters()) + 1; return (n > 0 && n <= 20); //fixme: totally arbitrary //fixme: should add condition that discrete-valued parents are fixed } } void FiniteSampler::burninOff() {}