Last active
March 3, 2023 07:57
-
-
Save xiangze/4c3c54d73bc7668ad139 to your computer and use it in GitHub Desktop.
Stick breaking process in stan
This file contains hidden or 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
#from BUGS book page 293 | |
#Stick-breaking process | |
#http://www2.imm.dtu.dk/courses/02443/projects/Roeder_JASA_1990.pdf | |
library(rstan) | |
C <- 10 | |
data <-read.csv("data.csv") | |
N <-nrow(data) | |
model <-stan_model("STB.stan") | |
fit <- sampling(model,v=data,N=N,C=C, iter = 1000, chains = 1) | |
traceplot(fit, ask=T) | |
print(fit, digit=2) |
This file contains hidden or 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
//from The BUGS book page 293 | |
//Stick-breaking process | |
data{ | |
int<lower=0> C;//num of cludter | |
int<lower=0> N;//data num | |
real v[N]; | |
} | |
parameters { | |
real amu; | |
real bmu; | |
vector[C] mu_mix; | |
vector[C] tau_mix; | |
real <lower=0,upper=1>q[C]; | |
} | |
transformed parameters{ | |
simplex [C] pi; | |
// simplex [C] p; | |
vector [C] mu_prec; | |
pi[1] <- q[1]; | |
for(j in 2:C) | |
pi[j] <- q[j]*(1-q[j-1])*pi[j-1]/q[j-1]; | |
for(j in 1:C){ | |
// pi[j] <- p[j]/sum(); | |
mu_prec[j] <- bmu*tau_mix[j]; | |
} | |
} | |
model { | |
// int<lower=0> group[C]; | |
real bprec; | |
real alpha; | |
real aprec; | |
amu ~ normal(0,0.001); | |
bmu ~ gamma(0.5,50); | |
bprec~gamma(2,1); | |
alpha <- 1; | |
aprec<-2; | |
for(j in 1:C){ | |
q[j]~beta(1,alpha); | |
mu_mix[j]~normal(amu,mu_prec[j]); | |
tau_mix[j]~gamma(aprec,bprec); | |
} | |
// for(i in 1:N){ | |
// group[i]~categorical(pi); | |
// for(j in 1:C) | |
// if(j==group[i]) | |
// gind[i,j]<-1; | |
// else | |
// gind[i,j]<-0; | |
// mu[i] <- mu.mix[group[i]]; | |
// tau[i] <- tau.mix[group[i]]; | |
// v[i]~normal(mu[i],tau[i]); | |
// } | |
for(i in 1:N){ | |
for(j in 1:C){ | |
increment_log_prob( | |
pi[j]+normal_log(v[i],mu_mix[j],tau_mix[j]) | |
); | |
} | |
} | |
} | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment