Skip to content

Commit 8e4f03d

Browse files
committed
Rel 0.1.0 - Switch to cmdstan 2.26.0
1 parent 435a08c commit 8e4f03d

12 files changed

Lines changed: 50 additions & 59 deletions

examples/arsenic/arsenic_logistic.stan

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,10 @@ transformed data {
1010
vector[p] mean_x;
1111
vector[p] sd_x;
1212
for (j in 1:p) {
13-
mean_x[j] <- mean(col(x,j));
14-
sd_x[j] <- sd(col(x,j));
13+
mean_x[j] = mean(col(x,j));
14+
sd_x[j] = sd(col(x,j));
1515
for (i in 1:N)
16-
z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
16+
z[i,j] = (x[i,j] - mean_x[j]) / sd_x[j];
1717
}
1818
}
1919

@@ -25,16 +25,17 @@ parameters {
2525

2626
model {
2727
vector[N] eta;
28-
eta <- beta0 + z*beta;
28+
beta0 ~ normal(0, phi);
2929
beta ~ normal(0, phi);
3030
phi ~ double_exponential(0, 10);
31+
eta = beta0 + z*beta;
3132
y ~ bernoulli_logit(eta);
3233
}
3334

3435
generated quantities {
3536
vector[N] log_lik;
3637
vector[N] eta;
37-
eta <- beta0 + z*beta;
38+
eta = beta0 + z*beta;
3839
for (i in 1:N)
39-
log_lik[i] <- bernoulli_logit_log(y[i],eta[i]);
40+
log_lik[i] = bernoulli_logit_lpmf(y[i] | eta[i]);
4041
}

examples/arsenic/arsenic_logistic_t.stan

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -13,12 +13,12 @@ transformed data {
1313
vector[p] mean_x;
1414
vector[p] sd_x;
1515
for (j in 1:p) {
16-
mean_x[j] <- mean(col(x,j));
17-
sd_x[j] <- sd(col(x,j));
16+
mean_x[j] = mean(col(x,j));
17+
sd_x[j] = sd(col(x,j));
1818
for (i in 1:N)
19-
z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
19+
z[i,j] = (x[i,j] - mean_x[j]) / sd_x[j];
2020
for (i in 1:Nt)
21-
zt[i,j] <- (xt[i,j] - mean_x[j]) / sd_x[j];
21+
zt[i,j] = (xt[i,j] - mean_x[j]) / sd_x[j];
2222
}
2323
}
2424
parameters {
@@ -28,7 +28,8 @@ parameters {
2828
}
2929
model {
3030
vector[N] eta;
31-
eta <- beta0 + z*beta;
31+
eta = beta0 + z*beta;
32+
beta0 ~ normal(0, phi);
3233
beta ~ normal(0, phi);
3334
phi ~ double_exponential(0, 10);
3435
y ~ bernoulli_logit(eta);
@@ -39,11 +40,11 @@ generated quantities {
3940
vector[Nt] log_likt;
4041
vector[N] eta;
4142
vector[Nt] etat;
42-
eta <- beta0 + z*beta;
43-
etat <- beta0 + zt*beta;
43+
eta = beta0 + z*beta;
44+
etat = beta0 + zt*beta;
4445
for (i in 1:N)
45-
log_lik[i] <- bernoulli_logit_log(y[i],eta[i]);
46+
log_lik[i] = bernoulli_logit_lpmf(y[i] | eta[i]);
4647
for (i in 1:Nt)
47-
log_likt[i] <- bernoulli_logit_log(yt[i],etat[i]);
48+
log_likt[i] = bernoulli_logit_lpmf(yt[i] | etat[i]);
4849
}
4950

examples/arsenic/compare.png

381 Bytes
Loading

examples/arsenic/demo_wells.jl

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,7 @@ using StatisticalRethinking
22
using JSON
33
using StanSample
44
using ParetoSmoothedImportanceSampling
5-
#using Statistics
65
using Printf
7-
#using StatsPlots
86

97
ProjDir = @__DIR__
108

@@ -110,6 +108,6 @@ end
110108

111109
# compare PSIS-LOO and k-fold-CV
112110
plot([-3.5, 0], [-3.5, 0], color=:red)
113-
scatter!(loos1, kfcvs, xlab = "PSIS-LOO", ylab = "10-fold-CV",
111+
scatter!(loos, kfcvs, xlab = "PSIS-LOO", ylab = "10-fold-CV",
114112
leg=false, color=:darkblue)
115113
savefig(joinpath(ProjDir, "compare.png"))

examples/roaches/diag_plot_1.png

-248 Bytes
Loading

examples/roaches/diag_plot_2.png

-552 Bytes
Loading

notebooks/arsenic.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -118,12 +118,12 @@ end
118118
md" ### k-fold-CV."
119119

120120
# ╔═╡ e86ee35c-5ebd-11eb-1932-ffbde638f11e
121-
md" ##### k-fold-CV should be used if several khats>0.5, in this case it is not needed, but provided as an example."
121+
md" ##### k-fold-CV should be used if several pk>0.5, in this case it is not needed, but provided as an example."
122122

123123
# ╔═╡ e86f0ab2-5ebd-11eb-1a0b-6958db3a7435
124124
begin
125125
model_str_2 = read(open(joinpath(ProjDir, "arsenic_logistic_t.stan")), String)
126-
sm3 = SampleModel("arsenic_logistic_t", model_str_2; tmpdir=tmpdir)
126+
sm3 = SampleModel("arsenic_logistic_t", model_str_2)
127127
end;
128128

129129
# ╔═╡ e70f73aa-5eb8-11eb-1960-bf6731681898

notebooks/cars.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
### A Pluto.jl notebook ###
2-
# v0.12.19
2+
# v0.12.20
33

44
using Markdown
55
using InteractiveUtils
@@ -23,7 +23,7 @@ md"
2323
# ╔═╡ 20d377b2-6008-11eb-364a-617b6934ecb2
2424
begin
2525
cd(psis_path)
26-
@quickactivate "PSIS"
26+
@quickactivate "ParetoSmoothedImportanceSampling"
2727
pkg"instantiate"
2828
end
2929

@@ -75,7 +75,7 @@ end
7575
if success(rc)
7676
nt_cars = read_samples(cars_stan_model);
7777
log_lik = nt_cars.log_lik'
78-
end
78+
end;
7979

8080
# ╔═╡ 20ed768a-6008-11eb-13f4-458ca1a29592
8181
begin

notebooks/roaches.jl

Lines changed: 8 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
### A Pluto.jl notebook ###
2-
# v0.12.18
2+
# v0.12.20
33

44
using Markdown
55
using InteractiveUtils
@@ -15,13 +15,13 @@ end
1515
# ╔═╡ d20c24f8-5ec2-11eb-3d45-d97fedebee8e
1616
begin
1717
cd(psis_path)
18-
@quickactivate "PSIS"
18+
@quickactivate "ParetoSmoothedImportanceSamplng"
1919
pkg"instantiate"
2020
end
2121

2222
# ╔═╡ e3552750-5e9f-11eb-324b-8df36d671c79
2323
begin
24-
ProjDir = joinpath(psis_path, "..", "Example", "roaches")
24+
ProjDir = joinpath(psis_path, "..", "examples", "roaches")
2525
df = CSV.read(joinpath(ProjDir, "roachdata.csv"), DataFrame)
2626
df.roach1 = df.roach1 / 100
2727
end;
@@ -93,13 +93,7 @@ if success(rc1)
9393
end
9494

9595
# ╔═╡ e3b82668-5e9f-11eb-1641-a9dfed9eb108
96-
begin
97-
scatter(pk1, xlab="Datapoint", ylab="Pareto shape k",
98-
marker=2.5, lab="Pk points")
99-
hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
100-
hline!([1], lab="pk = 1.0")
101-
title!("PSIS diagnostic plot for poisson-log model.")
102-
end
96+
pk_plot(pk1, title="PSIS diagnostic plot for poisson-log model.")
10397

10498
# ╔═╡ e3b8fc58-5e9f-11eb-0607-e5424a04df9c
10599
md" ##### Simple negative binomial regression example using the 2nd parametrization of the negative binomial distribution, see section 40.1-3 in the Stan reference guide."
@@ -165,21 +159,14 @@ end
165159

166160
# ╔═╡ ce509a0c-5ea8-11eb-3f2e-01023f29a1e3
167161
if success(rc2)
162+
168163
# Check the shape parameter k of the generalized Pareto distribution
169-
pk_good = sum(pk2 .<= 0.5)
170-
pk_ok = length(pk2[pk2 .<= 0.7]) - pk_good
171-
pk_bad = length(pk2[pk2 .<= 1]) - pk_good - pk_ok
172-
(good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk2 .> 1))
164+
165+
pk_qualify(pk2)
173166
end
174167

175168
# ╔═╡ e3f103f0-5e9f-11eb-159a-452961ed2619
176-
begin
177-
scatter(pk2, xlab="Datapoint", ylab="Pareto shape k",
178-
marker=2.5, lab="Pk points", leg=:topleft)
179-
hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
180-
hline!([1], lab="pk = 1.0")
181-
title!("PSIS diagnostic plot for neg-binomial model.")
182-
end
169+
pk_plot(pk2; title="PSIS diagnostic plot for neg-binomial model.", leg=:topright)
183170

184171
# ╔═╡ Cell order:
185172
# ╠═dcb4d418-5ec2-11eb-29d8-214f38b4d3ae

test/arsenic_logistic.stan

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,10 @@ transformed data {
1010
vector[p] mean_x;
1111
vector[p] sd_x;
1212
for (j in 1:p) {
13-
mean_x[j] <- mean(col(x,j));
14-
sd_x[j] <- sd(col(x,j));
13+
mean_x[j] = mean(col(x,j));
14+
sd_x[j] = sd(col(x,j));
1515
for (i in 1:N)
16-
z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
16+
z[i,j] = (x[i,j] - mean_x[j]) / sd_x[j];
1717
}
1818
}
1919

@@ -25,16 +25,17 @@ parameters {
2525

2626
model {
2727
vector[N] eta;
28-
eta <- beta0 + z*beta;
28+
beta0 ~ normal(0, phi);
2929
beta ~ normal(0, phi);
3030
phi ~ double_exponential(0, 10);
31+
eta = beta0 + z*beta;
3132
y ~ bernoulli_logit(eta);
3233
}
3334

3435
generated quantities {
3536
vector[N] log_lik;
3637
vector[N] eta;
37-
eta <- beta0 + z*beta;
38+
eta = beta0 + z*beta;
3839
for (i in 1:N)
39-
log_lik[i] <- bernoulli_logit_log(y[i],eta[i]);
40+
log_lik[i] = bernoulli_logit_lpmf(y[i] | eta[i]);
4041
}

0 commit comments

Comments
 (0)