voidGibbsSampling(int N, int K, VectorXi X, int MAXITER, int seed) {
// function to do Gibbs Sampling
// inputs:
// N: the number of data points
// K: the number of clusters
// X: the data
// MAXITER: the maximum number of iterations
// seed: the random seed value
// set random engine with the random seed value
std::default_random_engine engine(seed);
// set the output file
std::ofstream samples("GS-samples.csv");
// set the headers in the file
for (int k =0; k < K; k++) samples <<"a."<< k <<",";
for (int k =0; k < K; k++) samples <<"b."<< k <<",";
for (int k =0; k < K; k++) samples <<"lambda."<< k <<",";
for (int k =0; k < K; k++) samples <<"alpha."<< k <<",";
for (int k =0; k < K; k++) samples <<"pi."<< k <<",";
for (int n =0; n < N; n++) samples <<"s."<< n <<",";
samples <<"ELBO"<< std::endl;
// set variables
VectorXd a; // the shape parameter
VectorXd a_hat; // the modified shape parameter
VectorXd b; // the rate parameter
VectorXd b_hat; // the modified rate parameter
VectorXd alpha; // the concentration parameter
VectorXd alpha_hat; // the modified concentration parameter
VectorXd lambda; // the rate parameter
VectorXd pi; // the mixing parameter
MatrixXd eta(N,K); // the temporal parameter to sample s
VectorXi s; // the latent variable
MatrixXi S(N,K); // the latent variable (one-hot-labeling)
double ELBO; // ELBO
// set initial values
a = VectorXd::Constant(K,1,uniform_rng(0.1,2.0,engine));
b = VectorXd::Constant(K,1,uniform_rng(0.005,0.05,engine));
lambda = to_vector(gamma_rng(a,b,engine));
alpha = VectorXd::Constant(K,1,uniform_rng(10.0,200.0,engine));
pi = dirichlet_rng(alpha,engine);
// sampling
for (int i =1; i <= MAXITER; i++) {
// sample s and S
s = VectorXi::Zero(N,1); // initialize s with zeros
S = MatrixXi::Zero(N,K); // initialize S with zeros
for (int n =0; n < N; n++) {
eta.row(n) = X(n) * log(lambda) - lambda + log(pi);
eta.row(n) -= rep_row_vector(log_sum_exp(eta.row(n)),K);
s(n) = categorical_rng(exp(eta.row(n)), engine);
S(n,s(n)-1) =1;
}
// sample lambda
a_hat = a + S.transpose() * X;
b_hat = b + S.colwise().sum().transpose();
lambda = to_vector(gamma_rng(a_hat,b_hat,engine));
// sample pi
alpha_hat = alpha + S.colwise().sum().transpose();
pi = dirichlet_rng(alpha_hat,engine);
// calc ELBO
ELBO = calc_ELBO(N, K, X, a, b, alpha, a_hat, b_hat, alpha_hat);
// output
for (int k =0; k < K; k++) samples << a_hat(k) <<",";
for (int k =0; k < K; k++) samples << b_hat(k) <<",";
for (int k =0; k < K; k++) samples << lambda(k) <<",";
for (int k =0; k < K; k++) samples << alpha_hat(k) <<",";
for (int k =0; k < K; k++) samples << pi(k) <<",";
for (int n =0; n < N; n++) samples << s(n) <<",";
samples << ELBO << std::endl;
}
}
intmain(int argc, char*argv[]) {
// get inputs 1 ~ 4
std::string method = argv[1];
int N = atoi(argv[2]);
int K = atoi(argv[3]);
int seed = atoi(argv[4]);
if (method =="data") {
ただ、$\lambda_{1} = 44, \lambda_{2} = 77$ としてパラメータが渡されているデモデータと比較する上では、可視化の際などに色が反転してしまって見づらくなってしまうので、「# calculate EAP for lambda and sort them to control the order of clusters」の処理では、値を昇順に並び替えるように添字を振り直しています（対応してデモデータを生成するさいには常に昇順でパラメータを与えてやる必要があります）。
