Created
January 13, 2019 08:17
-
-
Save abikoushi/cf40ecb67c63421fce2e304cdb2363ba to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <RcppArmadillo.h> | |
using namespace Rcpp; | |
using namespace arma; | |
// [[Rcpp::depends(RcppArmadillo)]] | |
rowvec rdirichlet(const vec & alpha){ | |
int K = alpha.n_rows; | |
rowvec y(K); | |
for(int k=0;k<K ;k++){ | |
y[k] = R::rgamma(alpha[k],1); | |
} | |
rowvec out = y/sum(y); | |
return(out); | |
} | |
rowvec rcate(const rowvec & p){ | |
int K = p.n_cols; | |
rowvec cump = cumsum(p); | |
rowvec x(K); | |
x.fill(0); | |
double U = R::runif(0,1); | |
if(U<=cump[0]){ | |
x[0] = 1; | |
}else{ | |
for(int k=1; k<K; k++){ | |
if(cump[k-1]<U & U<=cump[k]){ | |
x[k] = 1; | |
} | |
} | |
} | |
return(x); | |
} | |
vec colSums(const mat & X){ | |
int nCols = X.n_cols; | |
vec out(nCols); | |
for(int i = 0; i < nCols; i++){ | |
out(i) = sum(X.col(i)); | |
} | |
return(out); | |
} | |
vec colMeans(const mat & X){ | |
int nCols = X.n_cols; | |
vec out(nCols); | |
for(int i = 0; i < nCols; i++){ | |
out(i) = sum(X.col(i))/X.n_rows; | |
} | |
return(out); | |
} | |
vec colSums(const imat & X){ | |
int nCols = X.n_cols; | |
vec out(nCols); | |
for(int i = 0; i < nCols; i++){ | |
out(i) = sum(X.col(i)); | |
} | |
return(out); | |
} | |
vec colMeans(const imat & X){ | |
int nCols = X.n_cols; | |
vec out(nCols); | |
for(int i = 0; i < nCols; i++){ | |
out(i) = sum(X.col(i))/X.n_rows; | |
} | |
return(out); | |
} | |
vec rowSums(const mat & X){ | |
int nRows = X.n_rows; | |
vec out(nRows); | |
for(int i = 0; i < nRows; i++){ | |
out(i) = sum(X.row(i)); | |
} | |
return(out); | |
} | |
rowvec softmax(const rowvec & x){ | |
double tmp = max(x); | |
rowvec res = exp(x-tmp)/sum(exp(x-tmp)); | |
return res; | |
} | |
double logsumexp(const rowvec & x){ | |
double tmp = max(x); | |
return tmp + std::log(sum(exp(x-tmp))); | |
} | |
double multi_lpmf(rowvec W, double M, rowvec p){ | |
int K = W.n_cols; | |
double lp = 0; | |
for(int k=0; k<K; k++){ | |
if(!(W[k]==0 && p[k]==0)){ | |
lp = lp + W[k]*std::log(p[k])-std::lgamma(W[k]+1); | |
} | |
} | |
lp = lp + std::lgamma(M+1); | |
return lp; | |
} | |
// [[Rcpp::export]] | |
List multmix_Gibbs(mat W, int L, vec phiini, mat pini, double alpha, int iter, double beta){ | |
vec M; | |
M = rowSums(W); | |
int N = W.n_rows; | |
int K = W.n_cols; | |
rowvec unnorm(L); | |
cube z = zeros(N,L,iter); | |
mat ztmp = zeros(N,L); | |
mat Wz = zeros(N,L); | |
mat phi = zeros(iter,L); | |
phi.row(0) = phiini.t(); | |
cube p = zeros(L,K,iter); | |
p.slice(0) = pini; | |
mat ptmp = pini; | |
mat loglik = zeros(iter-1); | |
for(int i=1;i<iter;i++){ | |
for(int n=0;n<N;n++){ | |
for(int l=0; l<L; l++){ | |
unnorm[l] = beta*(std::log(phi.row(i-1)[l]) + multi_lpmf(W.row(n),M[n],ptmp.row(l))); | |
} | |
loglik(i-1) += logsumexp(unnorm/beta); | |
ztmp.row(n) = rcate(softmax(unnorm)); | |
} | |
phi.row(i) = rdirichlet(beta*colSums(ztmp)+alpha); | |
for(int l=0; l<L; l++){ | |
Wz = W.each_col() % ztmp.col(l); | |
ptmp.row(l) = rdirichlet(beta*colSums(Wz)+alpha); | |
} | |
p.slice(i) = ptmp; | |
z.slice(i) = ztmp; | |
} | |
return List::create(Named("p")=p,_["phi"]=phi,_["z"]=z,_["loglik"]=loglik); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment