Commit 0f0d82fe authored by Jason Rhinelander's avatar Jason Rhinelander
Browse files

Did truncated t sampling as per R-Y et a. (2004)

This is about to get replaced, however, as it produces biased results!
parent 7382d976
......@@ -220,15 +220,16 @@ void Reader::interOptimize() {
auto max = demand_belief_.argmaxP(book->quality(), book->lifeSales(), book->age(), authored_books - 1, market_books, cost_unit);
const double &p = max.first;
const double &q = max.second;
#ifdef ERIS_DEBUG
if (p > 1000) {
ERIS_DBG("Found p=" << p << " > 1000");
ERIS_DBGVAR(id());
ERIS_DBGVAR(book->id());
ERIS_DBGVAR(p);
ERIS_DBGVAR(q);
ERIS_DBGVAR(demand_belief_.draw_discards);
ERIS_DBGVAR(demand_belief_.draw_success_cumulative);
ERIS_DBGVAR(demand_belief_.draw_discards_cumulative);
ERIS_DBGVAR(demand_belief_.draw_rejection_discards_last);
ERIS_DBGVAR(demand_belief_.draw_rejection_discards);
ERIS_DBGVAR(demand_belief_.draw_rejection_success);
std::cerr << "\n VectorXd beta(10); beta << ";
for (int i = 0; i < demand_belief_.beta().size(); i++) {
if (i > 0) std::cerr << ", ";
......@@ -250,6 +251,7 @@ void Reader::interOptimize() {
<< ", otherbooks = " << authored_books - 1 << ", marketbooks = " << market_books << "; } book;\n";
}
#endif
const double profit = (p - cost_unit) * q - cost_fixed;
if (profit > 0) {
......
#include "creativity/belief/Linear.hpp"
#include <eris/debug.hpp>
#include <eris/Random.hpp>
#include <Eigen/QR>
namespace creativity { namespace belief {
using namespace Eigen;
using eris::Random;
constexpr double Linear::NONINFORMATIVE_N, Linear::NONINFORMATIVE_S2, Linear::NONINFORMATIVE_Vc;
......@@ -13,11 +13,9 @@ Linear::Linear(
const Ref<const VectorXd> beta,
double s2,
const Ref<const MatrixXd> V,
double n,
std::shared_ptr<MatrixXd> V_inv,
std::shared_ptr<MatrixXd> V_chol_L
double n
)
: beta_(beta), s2_{s2}, V_(V), V_inv_{std::move(V_inv)}, V_chol_L_{std::move(V_chol_L)}, n_{n}, K_(beta_.rows())
: beta_(beta), s2_{s2}, V_(V.selfadjointView<Lower>()), n_{n}, K_(beta_.rows())
{
// Check that the given matrices conform
checkLogic();
......@@ -78,15 +76,21 @@ const MatrixXd& Linear::V() const { NO_EMPTY_MODEL; return V_; }
const MatrixXd& Linear::Vinv() const {
NO_EMPTY_MODEL;
if (not V_inv_)
V_inv_ = std::make_shared<MatrixXd>(V_.colPivHouseholderQr().inverse());
V_inv_ = std::make_shared<MatrixXd>(V_.colPivHouseholderQr().inverse().selfadjointView<Lower>());
return *V_inv_;
}
const MatrixXd& Linear::VcholL() const {
NO_EMPTY_MODEL;
if (not V_chol_L_)
V_chol_L_ = std::make_shared<MatrixXd>(V_.llt().matrixL());
V_chol_L_ = std::make_shared<MatrixXd>(V_.selfadjointView<Lower>().llt().matrixL());
return *V_chol_L_;
}
const MatrixXd& Linear::VcholLinv() const {
NO_EMPTY_MODEL;
if (not V_chol_L_inv_)
V_chol_L_inv_ = std::make_shared<MatrixXd>(VcholL().colPivHouseholderQr().inverse());
return *V_chol_L_inv_;
}
const bool& Linear::noninformative() const { NO_EMPTY_MODEL; return noninformative_; }
......@@ -102,26 +106,52 @@ const VectorXd& Linear::draw() {
if (last_draw_.size() != K_ + 1) last_draw_.resize(K_ + 1);
auto &rng = eris::Random::rng();
// beta is distributed as t(beta, s^2*V, n)
// That can be generated as beta + y*sqrt(n/q) where y ~ N(0, s^2*V), and q ~ chisq(n)
auto &rng = Random::rng();
// (beta,h) is distributed as a normal-gamma(beta, V, s2^{-1}, n), in Koop's Gamma distribution
// notation, or NG(beta, V, n/2, 2*s2^{-1}/n) in the more common G(shape,scale) notation
// (which std::gamma_distribution wants).
//
// Proof:
// Let $G_{k\theta}(k,\theta)$ be the shape ($k$), scale ($\theta$) notation. This has mean $k\theta$ and
// variance $k\theta^2$.
//
// Let $G_{Koop}(\mu,\nu)$ be Koop's notation, where $\mu$ is the mean and $\nu$ is the degrees of
// freedom, which has variance $\frac{2\mu^2}{\nu}$. Equating means and variances:
//
// \[
// k\theta = \mu
// k\theta^2 = \frac{2\mu^2}{\nu}
// \theta = \frac{2\mu}{\nu}
// k = \frac{2}{\nu}
// \]
// where the third equation follows from the first divided by the second, and fourth follows
// from the first divided by the third. Thus
// \[
// G_{Koop}(\mu,\nu) = G_{k\theta}(\frac{2}{\nu},\frac{2\mu}{\nu})
// \]
// To draw this, first draw a gamma-distributed "h" value (store its inverse)
last_draw_[K_] = 1.0 / std::gamma_distribution<double>(n_/2, 2/(s2_*n_))(rng);
// Now use that to draw a multivariate normal conditional on h, with mean beta and variance
// h^{-1} V; this is the beta portion of the draw:
last_draw_.head(K_) = multivariateNormal(beta_, VcholL(), std::sqrt(last_draw_[K_]));
// To generate y ~ N(0, SIGMA), generate N(0,1) and multiple by L from the cholesky
// decomposition of SIGMA, or in other words, s*L where LL'=V (so SIGMA=s^2*V)
VectorXd y(K_);
std::normal_distribution<double> stdnorm(0, 1);
for (unsigned int i = 0; i < K_; i++) y[i] = stdnorm(rng);
std::chi_squared_distribution<double> rchisqn(n_);
return last_draw_;
}
last_draw_.head(K_) = beta_ + sqrt(s2_ * n_ / rchisqn(rng)) * VcholL() * y;
VectorXd Linear::multivariateNormal(const Ref<const VectorXd> &mu, const Ref<const MatrixXd> &L, double s) {
if (mu.rows() != L.rows() or L.rows() != L.cols())
throw std::logic_error("multivariateNormal() called with non-conforming mu and L");
// h has distribution Gamma(n/2, 2/(n*s^2)), and s^2 is the inverse of h:
last_draw_[K_] = 1.0 / (std::gamma_distribution<double>(n_/2, 2/(s2_*n_))(rng));
// To draw such a normal, we need the lower-triangle Cholesky decomposition L of V, and a vector
// of K random \f$N(\mu=0, \sigma^2=h^{-1})\f$ values. Then \f$beta + Lz\f$ yields a \f$beta\f$
// draw of the desired distribution.
VectorXd z(mu.size());
for (unsigned int i = 0; i < z.size(); i++) z[i] = Random::rstdnorm();
return last_draw_;
return mu + L * (s * z);
}
const VectorXd& Linear::lastDraw() const {
......@@ -141,13 +171,23 @@ void Linear::discardForce(unsigned int burn) {
}
std::ostream& operator<<(std::ostream &os, const Linear &b) {
if (b.K() == 0)
os << "Linear model with no parameters (default constructed)";
else
os << "Linear model with " << b.K() << " parameters, beta_ =\n" << b.beta_;
b.print(os);
return os;
}
void Linear::print(std::ostream &os) const {
os << print_name();
if (K_ == 0) os << " model with no parameters (default constructed)";
else {
if (noninformative()) os << " (noninformative)";
os << " model: K=" << K_ << ", n=" << n_ << ", s2=" << s2_ <<
"\n beta = " << beta_.transpose().format(IOFormat(StreamPrecision, 0, ", ")) <<
"\n V = " << V_.format(IOFormat(6, 0, " ", "\n ")) << "\n";
}
}
std::string Linear::print_name() const { return "Linear"; }
void Linear::verifyParameters() const { NO_EMPTY_MODEL; }
// Called on an lvalue object, creates a new object with *this as prior
......@@ -194,10 +234,16 @@ void Linear::updateInPlace(const Ref<const VectorXd> &y, const Ref<const MatrixX
VectorXd beta_diff = beta_post - beta_;
beta_ = std::move(beta_post);
s2_ = (n_prior * s2_ + residualspost.squaredNorm() + beta_diff.transpose() * Vinv() * beta_diff) / n_;
//ERIS_DBG("orig nSSR: " << residualspost.squaredNorm()
/*ERIS_DBG("s2_ orig method: " << s2_);
ERIS_DBG("s2_ Koop: " << s2_alt);*/
V_inv_ = std::move(V_post_inv);
beta_ = std::move(beta_post);
if (V_chol_L_) V_chol_L_.reset(); // This will have to be recalculated
// The decompositions will have to be recalculated:
if (V_chol_L_) V_chol_L_.reset();
if (V_chol_L_inv_) V_chol_L_inv_.reset();
if (noninformative_) noninformative_ = false; // If we just updated a noninformative model, we aren't noninformative anymore
}
......@@ -231,14 +277,21 @@ void Linear::weakenInPlace(const double precision_scale) {
V_inv_ = std::make_shared<Eigen::MatrixXd>(*V_inv_ * precision_scale);
}
// Likewise for the Cholesky decomposition
// Likewise for the Cholesky decomposition (and its inverse)
if (V_chol_L_) {
if (V_chol_L_.unique())
*V_chol_L_ /= std::sqrt(precision_scale);
else
V_chol_L_ = std::make_shared<Eigen::MatrixXd>(*V_chol_L_ / std::sqrt(precision_scale));
}
if (V_chol_L_inv_) {
if (V_chol_L_inv_.unique())
*V_chol_L_inv_ *= std::sqrt(precision_scale);
else
V_chol_L_inv_ = std::make_shared<Eigen::MatrixXd>(*V_chol_L_inv_ * std::sqrt(precision_scale));
}
// And of course V gets scaled
V_ /= precision_scale;
return;
......
......@@ -8,9 +8,6 @@
namespace creativity { namespace belief {
/** Base class for a linear model with a natural conjugate, normal-gamma prior.
*
* FIXME: change (where appropriate) to use selfadjointView (esp. for Cholesky decomp, and possibly
* for inverse).
*/
class Linear {
public:
......@@ -50,33 +47,19 @@ class Linear {
*/
explicit Linear(unsigned int K);
// NB: if changing these constants, also change the above constructor documentation
static constexpr double
/** The value of `n` for a default noninformative model constructed using
* `Linear(unsigned int)`.
*/
//
NONINFORMATIVE_N = 1e-3,
/// The value of `s2` for a noninformative model constructed using `Linear(unsigned int)`
NONINFORMATIVE_S2 = 1.0,
/// The constant for the diagonals of the V matrix for a noninformative model
NONINFORMATIVE_Vc = 1e+8;
/** Constructs a Linear model with the given parameters. These parameters will be those
* used for the prior when updating.
*
* \param beta the coefficient mean parameters (which, because of restrictions, might not be
* the actual means).
*
* \param s2 the \f$\sigma^2\f$ value of the error term variance. Typically the \f$\sigma^2\f$ estimate.
*
* \param V the model's V matrix (where \f$s^2 V\f$ is the variance matrix of \f$\beta\f$).
* Note: only the lower triangle of the matrix will be used.
*
* \param n the number of data points supporting the other values (which can be a
* non-integer value).
* \param V_inv A shared pointer to the inverse of `V`, if already calculated. If the
* inverse has not already been calculated, it is better to omit this argument: the inverse
* will be calculated when needed.
* \param V_chol_L A shared pointer to the `L` matrix of the cholesky decomposition of `V`
* (where LL' = V). If the decomposition has not already been calculated, it is better to
* omit this argument: the decomposition will be calculated when needed.
*
* \throws std::runtime_error if any of (`K >= 1`, `V.rows() == V.cols()`, `K == V.rows()`)
* are not satisfied (where `K` is determined by the number of rows of `beta`).
......@@ -85,10 +68,8 @@ class Linear {
const Eigen::Ref<const Eigen::VectorXd> beta,
double s2,
const Eigen::Ref<const Eigen::MatrixXd> V,
double n,
std::shared_ptr<Eigen::MatrixXd> V_inv = nullptr,
std::shared_ptr<Eigen::MatrixXd> V_chol_L = nullptr
);
double n
);
/** Constructs a Linear model from std::vector<double>s containing the coefficients of beta
* and the lower triangle of V.
......@@ -125,6 +106,19 @@ class Linear {
/// Virtual destructor
virtual ~Linear() = default;
// NB: if changing these constants, also change the single-int, non-informative constructor documentation
static constexpr double
/** The value of `n` for a default noninformative model constructed using
* `Linear(unsigned int)`.
*/
//
NONINFORMATIVE_N = 1e-3,
/// The value of `s2` for a noninformative model constructed using `Linear(unsigned int)`
NONINFORMATIVE_S2 = 1.0,
/// The constant for the diagonals of the V matrix for a noninformative model
NONINFORMATIVE_Vc = 1e+8;
/** Virtual method called during construction to verify the model size. If this returns a
* non-zero value, the given parameters (beta, V for the regular constructor, K for the
* noninformative constructor) must agree with the returned value. If this returns 0, beta
......@@ -132,7 +126,6 @@ class Linear {
*/
virtual unsigned int fixedModelSize() const;
#define NO_EMPTY_MODEL if (K_ == 0) { throw std::logic_error("Cannot use default constructed model object as a model"); }
/** Accesses the base distribution means value of beta. Note that this is *not* necessarily
* the mean of beta and should not be used for prediction; rather it simply returns the
* distribution parameter value used, which may well not be the mean if any of the beta
......@@ -159,6 +152,10 @@ class Linear {
*/
const Eigen::MatrixXd& VcholL() const;
/** Accesses (calculating if not previous calculated) the inverse of `VcholL()`. Note that
* if VcholL() hasn't been calculated yet, this will calculate it. */
const Eigen::MatrixXd& VcholLinv() const;
/** Given a row vector of values \f$X^*\f$, predicts \f$y^*\f$ using the current model
* values. The default implementation provided by this class simply returns the mean \f$X^*
* \beta\f$ (the mean of the multivariate \f$t\f$ density for an unrestricted, natural
......@@ -172,10 +169,14 @@ class Linear {
*/
virtual double predict(const Eigen::Ref<const Eigen::RowVectorXd> &Xi);
/** Draws a vector of \f$\beta\f$ values and \f$s^2\f$ values distributed according to the
* model's structure. The default implementation simply returns the posterior beta and s^2
* values, which is suitable for a natural conjugate, normal-gamma distribution without
* any model restrictions.
/** Draws a vector of \f$\beta\f$ values and \f$h^{-1} = \sigma^2\f$ values distributed
* according to the model's parameters. The first `K()` values are the drawn \f$\beta\f$
* values, the last value is the drawn \f$h^{-1}\f$ value.
*
* In particular, this uses a gamma distribution to first draw an h value, then uses that h
* value to draw multivariate normal beta values. This means the \f$\beta\f$ values will have a
* multivariate t distribution with mean `beta()`, covariance parameter `s2()*V()`, and
* degrees of freedom parameter `n()`.
*
* \returns a const reference to the vector of values. This same vector is accessible by
* calling lastDraw(). Note that this vector is reused for subsequent draw() calls and so
......@@ -185,6 +186,36 @@ class Linear {
*/
virtual const Eigen::VectorXd& draw();
/** Draws a multivariate normal with mean \f$\mu\f$ covariance \f$s^2LL^\top\f$ (i.e. takes
* a constant and a Cholesky decomposition).
*
* \param mu the vector means
* \param L the Cholesky decomposition matrix
* \param s a standard deviation multiplier for the Cholesky decomposition matrix. Typically
* a \f$\sigma\f$ (NOT \f$\sigma^2\f$) value. If omitted, defaults to 1 (so that you can
* just pass the Cholesky decomposition of the full covariance matrix).
*
* \returns the random multivariate normal vector.
*
* \throws std::logic_error if mu and L have non-conforming sizes
*/
static Eigen::VectorXd multivariateNormal(
const Eigen::Ref<const Eigen::VectorXd> &mu,
const Eigen::Ref<const Eigen::MatrixXd> &L,
double s = 1.0);
/** Exception class thrown when draw() is unable to produce an admissable draw. Not thrown
* by this class (draws never fail) but available for subclass use.
*/
class draw_failure : public std::runtime_error {
public:
/** Constructor.
* \param what the exception message.
*/
draw_failure(const std::string &what) : std::runtime_error(what) {}
};
/** Returns a reference to the vector of \f$\beta\f$ and \f$s^2\f$ values generated by the
* last call to draw(). If draw() has not yet been called, the vector will be empty.
*/
......@@ -212,6 +243,16 @@ class Linear {
*/
friend std::ostream& operator << (std::ostream &os, const Linear &b);
/** Prints the Linear model to the given output stream. Called internally by operator<<,
* but subclassable. The model_base parameter is used for the first word of the output.
*/
virtual void print(std::ostream &os) const;
/** The display name of the model to use when printing it. Defaults to "Linear" but
* subclasses should override.
*/
virtual std::string print_name() const;
/** Using the calling object as a prior, uses the provided data to create a new Linear
* model.
*
......@@ -261,6 +302,7 @@ class Linear {
Linear weaken(double precision_scale) &&;
protected:
/** Weakens the current linear model. This functionality should only be used internally and
* by subclasses as required for move and copy update methods; weakening should be
* considered (externally) as a type of construction of a new object.
......@@ -305,11 +347,15 @@ class Linear {
*/
Eigen::MatrixXd V_;
/** The cached inverse of the prior V matrix, which isn't set until/unless needed. */
mutable std::shared_ptr<Eigen::MatrixXd> V_inv_;
mutable std::shared_ptr<Eigen::MatrixXd>
/// The cached inverse of the prior V matrix, which isn't set until/unless needed.
V_inv_,
/** The cached "L" matrix of the cholesky decomposition of V, where LL' = V. */
mutable std::shared_ptr<Eigen::MatrixXd> V_chol_L_;
/// The cached "L" matrix of the cholesky decomposition of V, where LL' = V.
V_chol_L_,
/// The cached inverse of the "L" matrix of the cholesky decomposition of V.
V_chol_L_inv_;
/// The number of data points supporting this model, which need not be an integer.
double n_;
......@@ -334,6 +380,4 @@ class Linear {
void checkLogic();
};
#undef NO_EMPTY_MODEL
}}
......@@ -2,10 +2,21 @@
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#include "creativity/belief/LinearRestricted.hpp"
#include <cmath>
#include <Eigen/QR>
#include <eris/algorithms.hpp>
#include <eris/Random.hpp>
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <eris/debug.hpp>
using namespace Eigen;
using eris::Random;
using boost::math::erfc_inv;
using boost::math::gamma_p;
using boost::math::gamma_q;
using boost::math::gamma_p_inv;
using boost::math::gamma_q_inv;
namespace creativity { namespace belief {
......@@ -47,6 +58,7 @@ void LinearRestricted::addRestriction(const Ref<const RowVectorXd> &R, double r)
restrict_values_[restrict_size_] = r;
restrict_select_.row(restrict_size_) = R;
restrict_size_++;
reset();
}
void LinearRestricted::addRestrictionGE(const Ref<const RowVectorXd> &R, double r) {
......@@ -61,6 +73,7 @@ void LinearRestricted::addRestrictions(const Ref<const MatrixXd> &R, const Ref<c
restrict_values_.segment(restrict_size_, num_restr) = r;
restrict_select_.middleRows(restrict_size_, num_restr) = R;
restrict_size_ += num_restr;
reset();
}
void LinearRestricted::addRestrictionsGE(const Ref<const MatrixXd> &R, const Ref<const VectorXd> &r) {
......@@ -69,6 +82,7 @@ void LinearRestricted::addRestrictionsGE(const Ref<const MatrixXd> &R, const Ref
void LinearRestricted::clearRestrictions() {
restrict_size_ = 0;
reset();
// Don't need to resize restrict_*, the values aren't used if restrict_size_ = 0
}
......@@ -85,53 +99,447 @@ void LinearRestricted::discardForce(unsigned int burn) {
void LinearRestricted::reset() {
Linear::reset();
mean_beta_draws_ = 0;
draw_discards = 0;
draw_discards_cumulative = 0;
draw_success_cumulative = 0;
draw_rejection_discards_last = 0;
draw_rejection_discards = 0;
draw_rejection_success = 0;
gibbs_D_.reset();
gibbs_alpha_.reset();
gibbs_last_.reset();
gibbs_draws_ = 0;
}
const Ref<const MatrixXd> LinearRestricted::R() const {
Ref<const MatrixXd> LinearRestricted::R() const {
return restrict_select_.topRows(restrict_size_);
}
const Ref<const VectorXd> LinearRestricted::r() const {
Ref<const VectorXd> LinearRestricted::r() const {
return restrict_values_.head(restrict_size_);
}
const VectorXd& LinearRestricted::draw() {
bool redraw = true;
draw_discards = 0;
while (redraw) {
return draw(draw_mode);
}
const VectorXd& LinearRestricted::draw(DrawMode m) {
// If they explicitly want rejection draw, do it:
if (m == DrawMode::Rejection)
return drawRejection();
// In auto mode we might try rejection sampling first
if (m == DrawMode::Auto) { // Need success rate < 0.1 to switch, with at least 20 samples.
long draw_rej_samples = draw_rejection_success + draw_rejection_discards;
double success_rate = draw_rejection_success / (double)draw_rej_samples;
if (success_rate < draw_auto_min_success_rate and draw_rej_samples >= draw_rejection_max_discards) {
// Too many failures, switch to Gibbs sampling
m = DrawMode::Gibbs;
}
else {
// Figure out how many sequential failed draws we'd need to hit the failure threshold,
// then try up to that many times.
long draw_tries = std::max<long>(
std::ceil(draw_rejection_success / draw_auto_min_success_rate),
draw_rejection_max_discards)
- draw_rej_samples;
try {
return drawRejection(draw_tries);
}
catch (draw_failure&) {
// Draw failure; switch to Gibbs
m = DrawMode::Gibbs;
}
}
}
// Either Gibbs was requested explicitly, or Auto was asked for and rejection sampling failed
// (either we tried it just now, or previous draws have too low a success rate), so use Gibbs
return drawGibbs();
}
void LinearRestricted::gibbsInitialize(const Ref<const VectorXd> &initial, unsigned long max_tries) {
constexpr double overshoot = 1.5;
if (initial.size() < K_ or initial.size() > K_+1) throw std::logic_error("LinearRestricted::gibbsInitialize() called with invalid initial vector (initial.size() != K())");
gibbs_draws_ = 0;
const MatrixXd &A = VcholLinv();
auto &rng = Random::rng();
if (restrict_size_ == 0) {
// No restrictions, nothing special to do!
if (not gibbs_last_ or gibbs_last_->size() != K_+1) gibbs_last_.reset(new VectorXd(K_+1));
gibbs_last_->head(K_) = A * initial.head(K_);
}
else {
// Otherwise we'll start at the initial value and update
VectorXd adjusted(initial.head(K_));
std::vector<size_t> violations;
violations.reserve(restrict_size_);
for (unsigned long trial = 1; ; trial++) {
violations.clear();
VectorXd v = restrict_select_.topRows(restrict_size_) * adjusted - restrict_values_.head(restrict_size_);
for (size_t i = 0; i < restrict_size_; i++) {
if (v[i] > 0) violations.push_back(i);
}
if (violations.size() == 0) break; // Restrictions satisfied!
else if (trial >= max_tries) { // Too many tries: give up
// Clear gibbs_last_ on failure (so that the next drawGibbs() won't try to use it)
gibbs_last_.reset();
throw constraint_failure("gibbsInitialize() couldn't find a way to satisfy the model constraints");
}
// Select a random constraint to fix:
size_t fix = violations[std::uniform_int_distribution<size_t>(0, violations.size()-1)(rng)];
// Aim at the nearest point on the boundary and overshoot (by 50%):
adjusted += -overshoot * v[fix] / restrict_select_.row(fix).squaredNorm() * restrict_select_.row(fix).transpose();
}
if (not gibbs_last_ or gibbs_last_->size() != K_+1) gibbs_last_.reset(new VectorXd(K_+1));
gibbs_last_->head(K_) = A * adjusted;
}
// Draw an appropriate sigma^2 value from the unconditional distribution, which should be
// suitable as an initial value (subsequent draws will be from the conditional distribution
// based on drawn beta values)
(*gibbs_last_)[K_] = 1.0 / std::gamma_distribution<double>(n_/2, 2/(s2_*n_))(rng);
}
VectorXd LinearRestricted::gibbsLast() {
return gibbs_last_ and gibbs_last_->size() == K_+1
? VcholL() * gibbs_last_->head(K_)
: VectorXd();
}
const VectorXd& LinearRestricted::drawGibbs() {
last_draw_mode = DrawMode::Gibbs;
const MatrixXd &A = VcholLinv();
const MatrixXd &Ainv = VcholL();
if (not gibbs_D_) gibbs_D_.reset(
restrict_size_ == 0
? new MatrixXd(0, K_) // It's rather pointless to use drawGibbs() with no restrictions, but allow it for debugging purposes
: new MatrixXd(restrict_select_.topRows(restrict_size_) * Ainv));
auto &D = *gibbs_D_;
int num_draws = 1;
if (gibbs_draws_ < draw_gibbs_burnin)
num_draws = draw_gibbs_burnin - gibbs_draws_ + 1;
else if (draw_gibbs_thinning > 1)
num_draws = draw_gibbs_thinning;
if (not gibbs_alpha_) gibbs_alpha_.reset(new VectorXd(A * beta_));
auto &alpha = *gibbs_alpha_;
auto &rng = Random::rng();
if (not gibbs_last_) {