182 lines
6.8 KiB
C
182 lines
6.8 KiB
C
#ifndef INCLUDE_GUARD_ISING_MCMC_H
|
|
#define INCLUDE_GUARD_ISING_MCMC_H
|
|
|
|
#include "R_api.h"
|
|
|
|
/** Sample a single value from the Ising model via a Monte-Carlo Markov-Chain
|
|
* method given the Ising model parameters in compact half-vectorized form.
|
|
*
|
|
* f(x) = p0(params) exp(vech(x x')' params)
|
|
*
|
|
* @important requires to be wrapped between `GetRNGstate()` and `PutRNGState()`
|
|
* calls. For example
|
|
*
|
|
* GetRNGstate()
|
|
* for (size_t sample = 0; sample < nr_samples; ++sample) {
|
|
* ising_mcmc_vech(warmup, dim, dim, params, &X[sample * dim]);
|
|
* }
|
|
* PutRNGState()
|
|
*/
|
|
static inline void ising_mcmc_vech(
|
|
const size_t warmup,
|
|
const size_t dim, const double* params,
|
|
int* X
|
|
) {
|
|
// Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
|
|
for (size_t i = 0, I = 0; i < dim; I += dim - i++) {
|
|
double invProb = 1.0 + exp(-params[I]);
|
|
X[i] = unif_rand() * invProb < 1.0;
|
|
}
|
|
|
|
// Skip the first samples of the Markov-Chain (warmup)
|
|
for (size_t skip = 0; skip < warmup; ++skip) {
|
|
// For every component
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
// Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
|
|
double log_odds = 0.0;
|
|
// Tracks position in half-vectorized storage
|
|
size_t J = i;
|
|
// Sum all log-odds for `j < i` (two way interactions)
|
|
for (size_t j = 0; j < i; J += dim - ++j) {
|
|
log_odds += X[j] ? params[J] : 0.0;
|
|
}
|
|
// Allways add log-odds for `i = j` (self-interaction)
|
|
log_odds += params[J];
|
|
// Continue adding log-odds for `j > i` (two way interactions)
|
|
J++;
|
|
for (size_t j = i + 1; j < dim; ++j, ++J) {
|
|
log_odds += X[j] ? params[J] : 0.0;
|
|
}
|
|
// Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
|
|
X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Thread save version of `ising_mcmc` (using thread save PRNG)
|
|
static inline void ising_mcmc_vech_thrd(
|
|
const size_t warmup,
|
|
const size_t dim, const double* params,
|
|
int* X, rng_seed_t seed
|
|
) {
|
|
// Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
|
|
for (size_t i = 0, I = 0; i < dim; I += dim - i++) {
|
|
double invProb = 1.0 + exp(-params[I]);
|
|
X[i] = unif_rand_thrd(seed) * invProb < 1.0;
|
|
}
|
|
|
|
// Skip the first samples of the Markov-Chain (warmup)
|
|
for (size_t skip = 0; skip < warmup; ++skip) {
|
|
// For every component
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
// Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
|
|
double log_odds = 0.0;
|
|
// Tracks position in half-vectorized storage
|
|
size_t J = i;
|
|
// Sum all log-odds for `j < i` (two way interactions)
|
|
for (size_t j = 0; j < i; J += dim - ++j) {
|
|
log_odds += X[j] ? params[J] : 0.0;
|
|
}
|
|
// Allways add log-odds for `i = j` (self-interaction)
|
|
log_odds += params[J];
|
|
// Continue adding log-odds for `j > i` (two way interactions)
|
|
J++;
|
|
for (size_t j = i + 1; j < dim; ++j, ++J) {
|
|
log_odds += X[j] ? params[J] : 0.0;
|
|
}
|
|
// Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
|
|
X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
/** Sample a single value from the Ising model via a Monte-Carlo Markov-Chain
|
|
* method given the Ising model parameters in symmetric matrix form `A`.
|
|
*
|
|
* f(x) = prob_0(A) exp(x' A x)
|
|
*
|
|
* The relation to the natural parameters `params` as used in the half-vectorized
|
|
* version `ising_mcmc_vech()`
|
|
*
|
|
* f(x) = prob_0(params) exp(vech(x x')' params)
|
|
*
|
|
* is illustrated via an example with `dim = 5`;
|
|
*
|
|
* p0 [ p0 p1/2 p2/2 p3/2 p4/2 ]
|
|
* p1 p5 [ p1/2 p5 p6/2 p7/2 p8/2 ]
|
|
* params = p2 p6 p9 => A = [ p2/2 p6/2 p9 p10/2 p11/2 ]
|
|
* p3 p7 p10 p12 [ p3/2 p7/2 p10/2 p12 p13/2 ]
|
|
* p4 p8 p11 p13 p14 [ p4/2 p8/2 p11/2 p13/2 p14 ]
|
|
*
|
|
* or the other way arroung given the symmetric matrix form `A` converted to
|
|
* the natural parameters
|
|
*
|
|
* [ a00 a10 a20 a30 a40 ] a00
|
|
* [ a10 a11 a21 a31 a41 ] 2*a10 a11
|
|
* A = [ a20 a21 a22 a32 a42 ] => params = 2*a20 2*a21 a22
|
|
* [ a30 a31 a32 a33 a43 ] 2*a30 2*a31 2*a32 a33
|
|
* [ a40 a41 a42 a43 a44 ] 2*a40 2*a41 2*a42 2*a43 a44
|
|
*
|
|
*/
|
|
static inline void ising_mcmc_mat(
|
|
const size_t warmup,
|
|
const size_t dim, const size_t ldA, const double* A,
|
|
int* X
|
|
) {
|
|
// Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
X[i] = unif_rand() * (1.0 + exp(-A[i * (ldA + 1)])) < 1.0;
|
|
}
|
|
|
|
// Skip the first samples of the Markov-Chain (warmup)
|
|
for (size_t skip = 0; skip < warmup; ++skip) {
|
|
// For every component
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
// Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
|
|
double log_odds = 0.0;
|
|
for (size_t j = 0; j < i; ++j) {
|
|
log_odds += (2 * X[j]) * A[i * ldA + j];
|
|
}
|
|
log_odds += A[i * (ldA + 1)];
|
|
for (size_t j = i + 1; j < dim; ++j) {
|
|
log_odds += (2 * X[j]) * A[i * ldA + j];
|
|
}
|
|
// Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
|
|
X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
// thread save version of `ising_mcmc_mat()` (using thread save PRNG)
|
|
static inline void ising_mcmc_mat_thrd(
|
|
const size_t warmup,
|
|
const size_t dim, const size_t ldA, const double* A,
|
|
int* X, rng_seed_t seed
|
|
) {
|
|
// Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
X[i] = unif_rand_thrd(seed) * (1.0 + exp(-A[i * (ldA + 1)])) < 1.0;
|
|
}
|
|
|
|
// Skip the first samples of the Markov-Chain (warmup)
|
|
for (size_t skip = 0; skip < warmup; ++skip) {
|
|
// For every component
|
|
for (size_t i = 0; i < dim; ++i) {
|
|
// Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
|
|
double log_odds = 0.0;
|
|
for (size_t j = 0; j < i; ++j) {
|
|
log_odds += (2 * X[j]) * A[i * ldA + j];
|
|
}
|
|
log_odds += A[i * (ldA + 1)];
|
|
for (size_t j = i + 1; j < dim; ++j) {
|
|
log_odds += (2 * X[j]) * A[i * ldA + j];
|
|
}
|
|
// Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
|
|
X[i] = unif_rand_thrd(seed) * (1.0 + exp(-log_odds)) < 1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
#endif /* INCLUDE_GUARD_ISING_MCMC_H */
|