Skip to content

Commit 9f273b4

Browse files
goedmangoedman
committed
Rel 0.0.1
1 parent 3cc6353 commit 9f273b4

12 files changed

Lines changed: 696 additions & 5 deletions

File tree

Example/arsenic_logistic.stan renamed to Example/arsenic/arsenic_logistic.stan

File renamed without changes.

Example/arsenic_logistic_t.stan renamed to Example/arsenic/arsenic_logistic_t.stan

File renamed without changes.

‎Example/arsenic/compare.png‎

19.9 KB
Loading

Example/cvit.jl renamed to Example/arsenic/cvit.jl

File renamed without changes.

Example/demo_wells.jl renamed to Example/arsenic/demo_wells.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -66,12 +66,12 @@ if success(rc2)
6666

6767
# Check the shape parameter k of the generalized Pareto distribution
6868
if all(pk .< 0.5)
69-
println("All Pareto k estimates OK (k < 0.5)")
69+
println("All Pareto k estimates OK (k < 0.5).")
7070
else
71-
pkn1 = sum((pk .>= 0.5) & (pk .< 1))
72-
pkn2 = sum(pk .>= 1)
73-
@printf(">> %d (%.0f%%) PSIS Pareto k estimates between 0.5 and 1\n", pkn1, pkn1/n*100)
74-
@printf(">> %d (%.0f%%) PSIS Pareto k estimates greater than 1\n", pkn2, pkn2/n*100)
71+
pk_good = sum(pk .<= 0.5)
72+
pk_ok = length(pk[pk .<= 0.7]) - pk_good
73+
pk_bad = length(pk[pk .<= 1]) - pk_good - pk_ok
74+
println((good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk .>= 1)))
7575
end
7676
end
7777

Example/wells.data.json renamed to Example/arsenic/wells.data.json

File renamed without changes.

‎Example/compare.png‎

-20.1 KB
Binary file not shown.

‎Example/roaches/roachdata.csv‎

Lines changed: 263 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,263 @@
1+
id,y,roach1,treatment,senior,exposure2
2+
1,153,308,1,0,0.8
3+
2,127,331.25,1,0,0.6
4+
3,7,1.67,1,0,1
5+
4,7,3,1,0,1
6+
5,0,2,1,0,1.14285714285714
7+
6,0,0,1,0,1
8+
7,73,70,1,0,0.8
9+
8,24,64.56,1,0,1.14285714285714
10+
9,2,1,0,0,1
11+
10,2,14,0,0,1.14285714285714
12+
11,0,138.25,0,0,1
13+
12,21,16,0,0,1
14+
13,0,97,0,0,1
15+
14,179,98,0,0,0.8
16+
15,136,44,0,0,1
17+
16,104,450,0,0,0.8
18+
17,2,36.67,0,0,0.8
19+
18,5,75,0,0,1
20+
19,1,2,0,0,1
21+
20,203,342.5,0,0,1
22+
21,32,22.5,0,0,1
23+
22,1,4,0,0,1
24+
23,135,94,0,0,1
25+
24,59,132.13,0,0,0.857142857142857
26+
25,29,80,0,0,1
27+
26,120,95,0,0,1
28+
27,44,34.13,0,0,0.8
29+
28,1,19,0,0,1
30+
29,2,25,0,0,1
31+
30,193,57.5,0,0,1
32+
31,13,3.11,0,0,1
33+
32,37,10.5,0,0,1
34+
33,2,47,0,0,1
35+
34,0,2.63,0,0,1
36+
35,3,204,0,0,1
37+
36,0,0,0,0,1
38+
37,0,0,0,0,1
39+
38,15,246.25,0,0,0.8
40+
39,11,76.22,0,0,1.57142857142857
41+
40,19,15,0,0,1.42857142857143
42+
41,0,14.88,0,0,1
43+
42,19,51,0,0,1
44+
43,4,15,0,0,1
45+
44,122,52.5,0,0,0.8
46+
45,48,13.75,0,0,1
47+
46,0,1.17,0,0,1
48+
47,0,1,0,0,1
49+
48,3,2,0,0,1
50+
49,0,0,0,0,1
51+
50,9,2,0,0,1
52+
51,0,2,0,0,1.14285714285714
53+
52,0,30,0,0,1.14285714285714
54+
53,0,7,0,0,1
55+
54,12,18,0,0,1
56+
55,0,2,0,0,0.8
57+
56,357,266,0,0,1
58+
57,11,174,1,0,0.914285714285714
59+
58,60,252,1,0,1
60+
59,0,0.88,1,0,1.28571428571429
61+
60,159,372.5,1,0,1
62+
61,50,42,1,0,1
63+
62,48,19,1,0,0.771428571428571
64+
63,178,263,1,0,1
65+
64,4,34,1,0,1.42857142857143
66+
65,6,1.75,1,0,1
67+
66,0,213.5,1,0,0.771428571428571
68+
67,33,13.13,1,0,1
69+
68,127,154,1,0,1
70+
69,4,21,1,0,1
71+
70,63,220,1,0,1
72+
71,88,38.32,1,0,1
73+
72,5,352.5,1,0,0.6
74+
73,0,7,1,0,1
75+
74,0,4,1,0,1
76+
75,62,59.5,1,0,1
77+
76,4,7.5,1,0,1
78+
77,150,112.88,1,0,1.28571428571429
79+
78,38,172,1,0,1
80+
79,0,13,1,0,1
81+
80,3,18,1,0,1
82+
81,1,0,1,0,1
83+
82,14,27,1,0,1
84+
83,77,148,1,0,1
85+
84,42,32,1,0,1.14285714285714
86+
85,21,28,1,0,1
87+
86,1,0,1,0,1
88+
87,45,28,1,0,1
89+
88,0,0,1,1,0.8
90+
89,0,14,1,1,0.8
91+
90,0,5,1,1,1
92+
91,0,0,1,1,1
93+
92,0,104,1,1,0.685714285714286
94+
93,183,27,1,1,1
95+
94,28,132,1,1,1
96+
95,49,258,1,1,1
97+
96,1,1,1,1,1
98+
97,0,2,1,1,1
99+
98,0,6,1,1,1
100+
99,3,3,1,1,0.8
101+
100,0,3,1,1,0.857142857142857
102+
101,0,0,1,1,1
103+
102,0,0,1,1,1
104+
103,0,0,1,1,1
105+
104,18,1.25,1,1,1
106+
105,0,0,1,1,1
107+
106,0,16,1,1,1
108+
107,5,68,1,1,0.4
109+
108,0,1,1,1,0.8
110+
109,19,18,1,1,1.14285714285714
111+
110,5,123.67,1,1,1
112+
111,0,2,1,1,0.8
113+
112,27,82.5,1,1,1
114+
113,0,0,1,1,0.2
115+
114,0,1.25,1,1,1
116+
115,77,171,1,1,1
117+
116,1,91.88,1,1,1
118+
117,3,5,1,1,1
119+
118,2,7,1,1,1
120+
119,0,0,1,1,1
121+
120,0,4,1,1,1
122+
121,22,53.75,1,1,1
123+
122,102,138.06,1,1,1.14285714285714
124+
123,0,1,1,1,1.14285714285714
125+
124,0,1.25,1,1,1
126+
125,0,28,1,1,0.8
127+
126,0,15,1,1,1
128+
127,0,0.88,1,1,0.8
129+
128,0,0,1,1,1
130+
129,0,33,1,1,1
131+
130,4,136.25,1,1,4.28571428571429
132+
131,12,127.5,1,1,0.8
133+
132,2,2,1,1,0.8
134+
133,0,2,1,1,0.8
135+
134,0,0,1,1,1
136+
135,1,3,1,1,1
137+
136,0,46,1,1,1.28571428571429
138+
137,40,68,1,1,1
139+
138,0,0,1,1,1
140+
139,1,49,1,1,1
141+
140,2,27.13,1,1,1
142+
141,27,45.5,1,0,1
143+
142,0,0,1,0,0.457142857142857
144+
143,2,4,1,0,1
145+
144,0,0,1,0,2.42857142857143
146+
145,0,1,1,0,1.42857142857143
147+
146,0,0,1,0,1
148+
147,0,0,1,0,1
149+
148,3,10.5,1,0,1.14285714285714
150+
149,1,0,1,0,1
151+
150,20,3,1,0,1.14285714285714
152+
151,0,0,1,0,1
153+
152,0,0,1,0,1
154+
153,0,0,1,0,1
155+
154,0,3.75,1,0,0.8
156+
155,0,0,1,0,1.14285714285714
157+
156,0,0,1,0,1
158+
157,0,0,0,0,0.8
159+
158,0,0.88,0,0,1.57142857142857
160+
159,0,2,1,0,1
161+
160,53,81,1,0,1.14285714285714
162+
161,69,31,1,0,2.28571428571429
163+
162,15,10,1,0,1
164+
163,0,10.23,1,0,0.571428571428571
165+
164,2,0,0,0,0.8
166+
165,4,13,1,0,0.857142857142857
167+
166,6,1,1,0,0.8
168+
167,8,33.25,1,0,0.857142857142857
169+
168,0,0,0,0,2.28571428571429
170+
169,0,53,0,0,0.8
171+
170,0,5,1,0,1
172+
171,18,157,1,0,0.857142857142857
173+
172,38,23.33,1,0,1
174+
173,0,6,1,0,1.02857142857143
175+
174,2,10,1,0,0.6
176+
175,18,100,1,0,0.857142857142857
177+
176,34,55,1,0,1
178+
177,1,0,1,0,0.8
179+
178,109,16.25,1,0,1
180+
179,5,7.78,1,0,1.48571428571429
181+
180,15,53,1,0,1
182+
181,0,2,1,0,1
183+
182,64,73,1,0,1
184+
183,0,0,1,0,0.8
185+
184,1,0,1,0,0.8
186+
185,0,0,1,0,1
187+
186,1,3.18,1,0,1
188+
187,3,5,0,0,0.857142857142857
189+
188,5,3,0,0,0.8
190+
189,7,0,0,0,1
191+
190,18,10,0,0,0.8
192+
191,1,1,0,0,1
193+
192,0,0,1,0,1
194+
193,0,1,1,0,0.8
195+
194,3,12,1,0,1
196+
195,3,17,1,0,1
197+
196,0,16,1,0,1
198+
197,19,2.5,1,0,1.42857142857143
199+
198,0,0,1,0,1.71428571428571
200+
199,8,21.88,1,0,0.771428571428571
201+
200,26,173,1,0,0.8
202+
201,50,111,1,0,1
203+
202,15,35,1,0,1.85714285714286
204+
203,0,0,1,0,1
205+
204,19,3,1,0,1
206+
205,5,2.1,1,0,1
207+
206,17,0,1,0,1
208+
207,121,49,1,0,1
209+
208,1,1.25,1,0,1
210+
209,0,1,1,0,0.8
211+
210,0,1,1,0,1
212+
211,0,0,1,0,1
213+
212,0,0,1,0,2.28571428571429
214+
213,4,54,1,0,0.8
215+
214,1,4.2,1,0,1.42857142857143
216+
215,14,51.25,1,0,1.02857142857143
217+
216,1,30,1,0,1.71428571428571
218+
217,25,196,0,0,1.14285714285714
219+
218,0,0,0,0,1.14285714285714
220+
219,14,2,0,0,1
221+
220,0,1.5,0,0,1
222+
221,59,96.25,0,0,1.14285714285714
223+
222,243,241.5,0,0,0.8
224+
223,80,140,0,0,1.42857142857143
225+
224,69,18,0,0,0.914285714285714
226+
225,14,3,0,0,1.14285714285714
227+
226,9,0.54,0,0,1.28571428571429
228+
227,38,82,0,0,0.6
229+
228,37,19,0,0,1
230+
229,48,18.75,0,0,0.8
231+
230,293,51,0,0,1.14285714285714
232+
231,7,0,0,0,1
233+
232,10,1,0,0,1.42857142857143
234+
233,19,0,0,0,1
235+
234,24,5.44,0,0,1
236+
235,91,0,0,0,1
237+
236,1,3,0,1,1
238+
237,0,0,0,1,1
239+
238,0,0,0,1,1
240+
239,0,0,0,1,1
241+
240,0,0,0,1,1
242+
241,148,28.75,0,1,1.14285714285714
243+
242,3,0,0,1,0.857142857142857
244+
243,26,3,0,1,1
245+
244,12,2,0,1,1
246+
245,77,135,0,1,1
247+
246,0,0,0,1,1
248+
247,7,68.25,0,1,1
249+
248,0,1,0,1,1
250+
249,1,0,0,1,1
251+
250,0,0,0,1,0.6
252+
251,17,0,0,1,1.02857142857143
253+
252,0,1,0,1,1
254+
253,7,1,0,1,1
255+
254,11,2.5,0,1,0.685714285714286
256+
255,6,51.25,0,1,0.8
257+
256,50,13.13,0,1,1
258+
257,1,0,0,1,0.8
259+
258,0,0,0,1,1.48571428571429
260+
259,0,0,0,1,1
261+
260,0,0,0,1,1
262+
261,171,0,0,1,1
263+
262,8,0,0,1,1

‎Example/roaches/roaches.jl‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
using StatisticalRethinking, StanSample, PSIS
2+
using Printf
3+
4+
ProjDir = @__DIR__
5+
6+
df = CSV.read(joinpath(ProjDir, "roachdata.csv"), DataFrame)
7+
df.roach1 = df.roach1 / 100
8+
9+
roaches1_stan = "
10+
data {
11+
int<lower=0> N;
12+
vector[N] exposure2;
13+
vector[N] roach1;
14+
vector[N] senior;
15+
vector[N] treatment;
16+
int y[N];
17+
}
18+
transformed data {
19+
vector[N] log_expo;
20+
21+
log_expo = log(exposure2);
22+
}
23+
parameters {
24+
vector[4] beta;
25+
}
26+
model {
27+
y ~ poisson_log(log_expo + beta[1] + beta[2] * roach1 + beta[3] * treatment
28+
+ beta[4] * senior);
29+
}
30+
generated quantities {
31+
vector[N] log_lik;
32+
vector[N] eta;
33+
eta = log_expo + beta[1] + beta[2] * roach1 + beta[3] * treatment
34+
+ beta[4] * senior;
35+
for (i in 1:N)
36+
log_lik[i] <- poisson_log_lpmf(y[i] | eta[i]);
37+
}
38+
";
39+
40+
data = (N = size(df, 1), y = Int.(df.y), roach1=df.roach1,
41+
exposure2=df.exposure2, senior=df.senior, treatment=df.treatment)
42+
n = size(df, 1)
43+
44+
tmpdir=joinpath(ProjDir, "tmp")
45+
sm1 = SampleModel("roaches1", roaches1_stan; tmpdir)
46+
rc1 = stan_sample(sm1; data)
47+
48+
if success(rc1)
49+
#roaches1_df = read_samples(sm1; output_format=:dataframe)
50+
#precis(roaches1_df)
51+
nt1 = read_samples(sm1)
52+
53+
# Compute LOO and standard error
54+
log_lik1 = nt1.log_lik'
55+
loo1, loos1, pk1 = psisloo(log_lik1)
56+
elpd_loo1 = sum(loos1)
57+
se_elpd_loo1 = std(loos1) * sqrt(n)
58+
@printf(">> elpd_loo = %.1f, SE(elpd_loo) = %.1f\n", elpd_loo1, se_elpd_loo1)
59+
60+
# Check the shape parameter k of the generalized Pareto distribution
61+
if all(pk1 .< 0.5)
62+
println("All Pareto k estimates OK (k < 0.5).")
63+
else
64+
pk_good = sum(pk1 .<= 0.5)
65+
pk_ok = length(pk1[pk1 .<= 0.7]) - pk_good
66+
pk_bad = length(pk1[pk1 .<= 1]) - pk_good - pk_ok
67+
println((good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk1 .> 1)))
68+
end
69+
end
70+
71+
begin
72+
scatter(pk1, xlab="Datapoint", ylab="Pareto shape k",
73+
marker=2.5, lab="Pk points")
74+
hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
75+
hline!([1], lab="pk = 1.0")
76+
title!("PSIS diagnostic plot")
77+
end

0 commit comments

Comments
 (0)