Skip to content

Commit 1b95e9e

Browse files
sycliksyclik
committed
Merge pull request #154 from stan-dev/feature/update-bugs-examples
Feature/update bugs examples: vectorization, use bernoulli/binomial_logit
2 parents ad2c551 + fe8d24f commit 1b95e9e

40 files changed

+284
-393
lines changed

‎bugs_examples/vol1/dogs/makefile‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ LIBFLAGS = -L$(STAN_HOME)/bin -lstan
99
$(PGM) :
1010
$(STAN_HOME)/bin/stanc --name=$(PGM) $(PGM).stan
1111
$(CXX) -Wall -O3 -DNDEBUG $(CPPFLAGS) $(PGM).cpp -o $(PGM) $(LIBFLAGS)
12-
./$(PGM) --data=$(PGM).data.R --seed=1340384924 --chain_id=4 # --iters=20000
12+
./$(PGM) --data=$(PGM).data.R --seed=1340384924 --chain_id=4 --diagnostics=d.csv# --iters=20000
1313

1414
clean :
15-
rm $(PGM).cpp samples.csv $(PGM)
15+
rm -rf $(PGM).cpp samples.csv $(PGM)

‎bugs_examples/vol1/dogs/post.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11

22
library(coda)
3-
post <- read.csv(file = "samples.csv", header = TRUE, comment.char = '#')[, -(1:3)]
3+
post <- read.csv(file = "samples.csv", header = TRUE, comment.char = '#')[, -(1:4)]
44
summary(as.mcmc(post))
55
plot(as.mcmc(post))
66

‎bugs_examples/vol1/epil/epil.stan‎

Lines changed: 27 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -10,16 +10,31 @@ data {
1010
int<lower=0> y[N, T];
1111
int<lower=0> Trt[N];
1212
int<lower=0> V4[T];
13-
real log_Base4[N];
14-
real log_Age[N];
15-
real BT[N];
13+
vector[N] log_Base4;
14+
vector[N] log_Age;
15+
vector[N] BT;
1616
real log_Age_bar;
1717
real Trt_bar;
1818
real BT_bar;
1919
real V4_bar;
2020
real log_Base4_bar;
2121
}
2222

23+
transformed data {
24+
vector[T] V4_c;
25+
vector[N] log_Base4_c;
26+
vector[N] log_Age_c;
27+
vector[N] BT_c;
28+
vector[N] Trt_c;
29+
log_Base4_c <- log_Base4 - log_Base4_bar;
30+
log_Age_c <- log_Age - log_Age_bar;
31+
BT_c <- BT - BT_bar;
32+
for (i in 1:T)
33+
V4_c[i] <- V4[i] - V4_bar;
34+
for (i in 1:N)
35+
Trt_c[i] <- Trt[i] - Trt_bar;
36+
}
37+
2338
parameters {
2439
real a0;
2540
real alpha_Base;
@@ -28,7 +43,7 @@ parameters {
2843
real alpha_Age;
2944
real alpha_V4;
3045
real b1[N];
31-
real b[N, T];
46+
vector[T] b[N];
3247
real<lower=0> sigmasq_b;
3348
real<lower=0> sigmasq_b1;
3449
}
@@ -40,6 +55,7 @@ transformed parameters {
4055
sigma_b1 <- sqrt(sigmasq_b1);
4156
}
4257

58+
4359
model {
4460
a0 ~ normal(0, 100);
4561
alpha_Base ~ normal(0, 100);
@@ -49,17 +65,14 @@ model {
4965
alpha_V4 ~ normal(0, 100);
5066
sigmasq_b1 ~ inv_gamma(.001, .001);
5167
sigmasq_b ~ inv_gamma(.001, .001);
68+
b1 ~ normal(0, sigma_b1);
5269
for(n in 1:N) {
53-
b1[n] ~ normal(0, sigma_b1);
54-
for(t in 1:T) {
55-
b[n, t] ~ normal(0, sigma_b);
56-
y[n, t] ~ poisson(exp(a0 + alpha_Base * (log_Base4[n] - log_Base4_bar)
57-
+ alpha_Trt * (Trt[n] - Trt_bar)
58-
+ alpha_BT * (BT[n] - BT_bar)
59-
+ alpha_Age * (log_Age[n] - log_Age_bar)
60-
+ alpha_V4 * (V4[t] - V4_bar)
61-
+ b1[n] + b[n, t]));
62-
}
70+
b[n] ~ normal(0, sigma_b);
71+
y[n] ~ poisson_log(a0 + alpha_Base * log_Base4_c[n]
72+
+ alpha_Trt * Trt_c[n]
73+
+ alpha_BT * BT_c[n]
74+
+ alpha_Age * log_Age_c[n]
75+
+ b1[n] + alpha_V4 * V4_c + b[n]);
6376
}
6477
}
6578

‎bugs_examples/vol1/equiv/equiv.stan‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -33,10 +33,11 @@ transformed parameters {
3333
}
3434

3535
model {
36-
for (p in 1:P) {
37-
for (n in 1:N) {
38-
Y[n, p] ~ normal(mu + sign[T[n, p]] * phi / 2 + sign[p] * pi / 2 + delta[n], sigma1);
39-
}
36+
for (n in 1:N) {
37+
vector[P] m;
38+
for (p in 1:P)
39+
m[p] <- mu + sign[T[n, p]] * phi / 2 + sign[p] * pi / 2 + delta[n];
40+
Y[n] ~ normal(m, sigma1);
4041
}
4142
delta ~ normal(0, sigma2);
4243
sigmasq1 ~ inv_gamma(.001, .001);

‎bugs_examples/vol1/leuk/leuk.stan‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,6 @@ model {
3535
for(j in 1:NT) {
3636
dL0[j] ~ gamma(r * (t[j + 1] - t[j]) * c, c);
3737
for(i in 1:N) {
38-
// lp__ <- lp__ + if_else(Y[i, j], poisson_log(dN[i, j], Y[i, j] * exp(beta * Z[i]) * dL0[j]), 0);
3938
if (Y[i, j] != 0) lp__ <- lp__ + poisson_log(dN[i, j], Y[i, j] * exp(beta * Z[i]) * dL0[j]);
4039
}
4140
}

‎bugs_examples/vol1/magnesium/magnesium.stan‎

Lines changed: 9 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ transformed data {
2626
}
2727

2828
parameters {
29-
real mu[N_priors];
29+
real<lower=-10,upper=10> mu[N_priors];
3030
real theta[N_priors, N_studies];
3131
real<lower=0,upper=1> pc[N_priors, N_studies];
3232
real<lower=0> inv_tau_sq_1;
@@ -67,23 +67,19 @@ model {
6767
// Prior 6: Half-Normal on tau.sqrd
6868
tau_sq_6 ~ normal(0, p0_sigma) T[0,];
6969

70-
for (prior in 1:N_priors)
71-
mu[prior] ~ uniform(-10, 10);
70+
mu ~ uniform(-10, 10);
7271

7372
for (prior in 1:N_priors) {
74-
for (study in 1:N_studies) {
75-
pc[prior, study] ~ uniform(0,1);
76-
theta[prior, study] ~ normal(mu[prior], tau[prior]);
77-
}
73+
pc[prior] ~ uniform(0,1);
74+
theta[prior] ~ normal(mu[prior], tau[prior]);
7875
}
7976
for (prior in 1:N_priors) {
80-
for (study in 1:N_studies) {
81-
rc[study] ~ binomial(nc[study], pc[prior, study]);
82-
rt[study] ~ binomial(nt[study], inv_logit(theta[prior, study] +
83-
logit(pc[prior, study])));
84-
}
77+
vector[N_studies] tmpm;
78+
for (i in 1:N_studies)
79+
tmpm[i] <- theta[prior, i] + logit(pc[prior, i]);
80+
rc ~ binomial(nc, pc[prior]);
81+
rt ~ binomial_logit(nt, tmpm);
8582
}
86-
8783
}
8884
generated quantities {
8985
real OR[N_priors];

‎bugs_examples/vol1/oxford/oxford.stan‎

Lines changed: 9 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -7,33 +7,29 @@ data {
77
int<lower=0> n1[K];
88
int<lower=0> r0[K];
99
int<lower=0> r1[K];
10-
int year[K];
10+
vector[K] year;
1111
}
1212
transformed data {
13-
int yearsq[K];
14-
for (i in 1:K)
15-
yearsq[i] <- year[i] * year[i];
13+
vector[K] yearsq;
14+
yearsq <- year .* year;
1615
}
1716
parameters {
18-
real mu[K];
17+
vector[K] mu;
1918
real alpha;
2019
real beta1;
2120
real beta2;
2221
real<lower=0> sigma_sq;
23-
real b[K];
22+
vector[K] b;
2423
}
2524
transformed parameters {
2625
real<lower=0> sigma;
2726
sigma <- sqrt(sigma_sq);
2827
}
2928
model {
30-
for (i in 1:K) {
31-
r0[i] ~ binomial(n0[i], inv_logit(mu[i]));
32-
r1[i] ~ binomial(n1[i],
33-
inv_logit(mu[i] + alpha + beta1 * year[i] + beta2 * (yearsq[i] - 22) + sigma * b[i]));
34-
b[i] ~ normal(0, 1);
35-
mu[i] ~ normal(0, 1000);
36-
}
29+
r0 ~ binomial_logit(n0, mu);
30+
r1 ~ binomial_logit(n1, alpha + mu + beta1 * year + beta2 * (yearsq - 22) + b * sigma);
31+
b ~ normal(0, 1);
32+
mu ~ normal(0, 1000);
3733

3834
alpha ~ normal(0.0, 1000);
3935
beta1 ~ normal(0.0, 1000);

‎bugs_examples/vol1/salm/salm.data.R‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
1-
doses <-
1+
Ndoses <-
22
6
3-
plates <-
3+
Nplates <-
44
3
55
y <-
66
structure(c(15, 16, 16, 27, 33, 20, 21, 18, 26, 41, 38, 27, 29,

‎bugs_examples/vol1/salm/salm.init.R‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,8 @@ gamma <- 0
44
tau <- 0.1
55
lambda <-
66
structure(c(0, 0, 0,
7-
0, 0, 0,
8-
0, 0, 0,
9-
0, 0, 0,
10-
0, 0, 0,
11-
0, 0, 0),
12-
.Dim = c(6, 3))
7+
0, 0, 0,
8+
0, 0, 0,
9+
0, 0, 0,
10+
0, 0, 0,
11+
0, 0, 0), .Dim = c(6, 3))

‎bugs_examples/vol1/salm/salm.stan‎

Lines changed: 18 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1,34 +1,34 @@
11
## http://www.openbugs.info/Examples/Salm.html
22
## this matches the jags implementation
33
data {
4-
int<lower=0> doses;
5-
int<lower=0> plates;
6-
int<lower=0> y[doses,plates];
7-
real x[doses];
4+
int<lower=0> Ndoses;
5+
int<lower=0> Nplates;
6+
int<lower=0> y[Ndoses,Nplates];
7+
real x[Ndoses];
88
}
99
transformed data {
10-
real logx[doses];
10+
real logx[Ndoses];
1111
real mean_x;
1212
real mean_logx;
13-
real centered_x[doses];
14-
real centered_logx[doses];
13+
real centered_x[Ndoses];
14+
real centered_logx[Ndoses];
1515

1616
mean_x <- mean(x);
17-
for (dose in 1:doses)
17+
for (dose in 1:Ndoses)
1818
centered_x[dose] <- x[dose] - mean_x;
1919

20-
for (dose in 1:doses)
20+
for (dose in 1:Ndoses)
2121
logx[dose] <- log(x[dose] + 10);
2222
mean_logx <- mean(logx);
23-
for (dose in 1:doses)
23+
for (dose in 1:Ndoses)
2424
centered_logx[dose] <- logx[dose] - mean_logx;
2525
}
2626
parameters {
2727
real alpha_star;
2828
real beta;
2929
real gamma;
3030
real<lower=0> tau;
31-
real lambda[doses,plates];
31+
vector[Nplates] lambda[Ndoses];
3232
}
3333
transformed parameters {
3434
real<lower=0> sigma;
@@ -37,18 +37,17 @@ transformed parameters {
3737
alpha <- alpha_star - beta * mean_logx - gamma * mean_x;
3838
sigma <- 1.0 / sqrt(tau);
3939
}
40+
4041
model {
4142
alpha_star ~ normal(0.0,1.0E3);
4243
beta ~ normal(0.0,1000);
4344
gamma ~ normal(0.0,1000);
4445
tau ~ gamma(0.001,0.001);
45-
for (dose in 1:doses) {
46-
for (plate in 1:plates) {
47-
lambda[dose,plate] ~ normal(0.0, sigma);
48-
y[dose,plate] ~ poisson(exp(alpha_star
49-
+ beta * centered_logx[dose]
50-
+ gamma * centered_x[dose]
51-
+ lambda[dose,plate]) );
52-
}
46+
for (dose in 1:Ndoses) {
47+
lambda[dose] ~ normal(0.0, sigma);
48+
y[dose] ~ poisson_log(alpha_star
49+
+ beta * centered_logx[dose]
50+
+ gamma * centered_x[dose]
51+
+ lambda[dose]);
5352
}
5453
}

0 commit comments

Comments
 (0)