Rel 7.0.0

Author: goedman

Examples/Basic_intros/chimpanzees.jl

@@ -0,0 +1,119 @@
+# A walk-through example (using StanSample.jl)
+
+using StanSample, Distributions
+
+N = 100
+df = DataFrame(
+    height = rand(Normal(10, 2), N),
+    leg_prop = rand(Uniform(0.4, 0.5), N),
+)
+df.leg_left = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
+df.leg_right = df.leg_prop .* df.height + rand(Normal(0, 0.03), N)
+
+stan6_1 = "
+data {
+  int <lower=1> N;
+  vector[N] H;
+  vector[N] LL;
+  vector[N] LR;
+}
+parameters {
+  real a;
+  vector[2] b;
+  real <lower=0> sigma;
+}
+model {
+  vector[N] mu;
+  mu = a + b[1] * LL + b[2] * LR;
+  a ~ normal(10, 100);
+  b ~ normal(2, 10);
+  sigma ~ exponential(1);
+  H ~ normal(mu, sigma);
+}
+";
+
+m6_1s = SampleModel("m6.1s", stan6_1);
+data = (H = df.height, LL = df.leg_left, LR = df.leg_right, N = size(df, 1))
+rc6_1s = stan_sample(m6_1s; data);
+
+if success(rc6_1s)
+    chns = read_samples(m6_1s)
+
+    # Display the chns
+
+    chns |> display
+    println()
+
+    # Display the keys
+
+    axiskeys(chns) |> display
+    println()
+
+    chns(chain=1) |> display
+    println()
+
+    chns[:, 1, 1] |> display
+    println()
+
+    chns(chain=1, param=:a) |> display
+    println()
+
+    chns(chain=[1, 3], param=[:a, :sigma]) |> display
+    println()
+
+    # Select all elements starting with 'a'
+
+    chns_b = matrix(chns, :b)
+    chns_b |> display
+    println()
+
+    mean(chns_b, dims=1) |> display
+    println()
+
+    typeof(chns_b.data) |> display
+    println()
+
+    ndraws_b, nchains_b, nparams_b = size(chns_b)
+    chn_b = reshape(chns_b, ndraws_b*nchains_b, nparams_b)
+    println()
+
+    for row in eachrow(chn_b)
+        # ...
+    end
+
+    # Or use read_samples to only use chains 2 and 4 using the chains kwarg.
+
+    chns2 = read_samples(m6_1s; chains=[2, 4])
+    chns2_b = matrix(chns2, :b)
+    ndraws2_b, nchains2_b, nparams2_b = size(chns2_b)
+    chn2_b = reshape(chns2_b, ndraws2_b*nchains2_b, nparams2_b)
+    mean(chns2_b, dims=1) |> display
+
+end
+
+# Some more options:
+
+init = (a = 2.0, b = [1.0, 2.0], sigma = 1.0)
+rc6_2s = stan_sample(m6_1s; data, init);
+
+if success(rc6_2s)
+
+    read_summary(m6_1s, true)
+    println()
+
+    post6_1s_df = read_samples(m6_1s, :dataframe)
+    post6_1s_df |> display
+    println()
+
+    part6_1s = read_samples(m6_1s, :particles)
+    part6_1s |> display
+    println()
+
+    nt6_1s = read_samples(m6_1s, :namedtuple)
+    nt6_1s.b |> display
+    println()
+
+    mean(nt6_1s.b, dims=2) |> display
+    println()
+
+end

Examples/Basic_intros/chimpanzees_keyedarray.jl

@@ -0,0 +1,100 @@
+# Load Julia packages (libraries)
+
+using Pkg, DrWatson
+
+using StanSample
+
+df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), DataFrame);
+
+# Define the Stan language model
+
+stan10_4 = "
+data{
+    int N;
+    int N_actors;
+    int pulled_left[N];
+    int prosoc_left[N];
+    int condition[N];
+    int actor[N];
+}
+parameters{
+    vector[N_actors] a;
+    real bp;
+    real bpC;
+}
+model{
+    vector[N] p;
+    bpC ~ normal( 0 , 10 );
+    bp ~ normal( 0 , 10 );
+    a ~ normal( 0 , 10 );
+    for ( i in 1:504 ) {
+        p[i] = a[actor[i]] + (bp + bpC * condition[i]) * prosoc_left[i];
+        p[i] = inv_logit(p[i]);
+    }
+    pulled_left ~ binomial( 1 , p );
+}
+";
+
+data = (N = size(df, 1), N_actors = length(unique(df.actor)), 
+    actor = df.actor, pulled_left = df.pulled_left,
+    prosoc_left = df.prosoc_left, condition = df.condition);
+
+# Sample using cmdstan
+
+m10_4s = SampleModel("m10.4s", stan10_4)
+rc10_4s = stan_sample(m10_4s; data);
+
+if success(rc10_4s)
+    chns = read_samples(m10_4s)
+
+    # Display the chns
+
+    chns |> display
+    println()
+
+    # Display the keys
+
+    axiskeys(chns) |> display
+    println()
+
+    mean(chns_a, dims=1) |> display
+    println()
+
+    chns(chain=1) |> display
+    println()
+
+    chns[:, 1, 8] |> display
+    println()
+
+    chns(chain=1, param=:bp) |> display
+    println()
+
+    chns(chain=[1, 3], param=[:bp, :bpC]) |> display
+    println()
+
+    # Select all elements starting with 'a'
+
+    chns_a = matrix(chns, :a)
+    chns_a |> display
+    println()
+
+    typeof(chns_a.data) |> display
+    println()
+
+    ndraws_a, nchains_a, nparams_a = size(chns_a)
+    chn_a = reshape(chns_a, ndraws_a*nchains_a, nparams_a)
+    println()
+
+    for row in eachrow(chn_a)
+        # ...
+    end
+
+    # Or use read_samples to only use chains 2 and 4 using the chains kwarg.
+
+    chns2 = read_samples(m10_4s; chains=[2, 4])
+    chns2_a = matrix(chns2, :a)
+    ndraws2_a, nchains2_a, nparams2_a = size(chns2_a)
+    chn2_a = reshape(chns2_a, ndraws2_a*nchains2_a, nparams2_a)
+    mean(chns2_a, dims=1) |> display
+
+end

Examples/Basic_intros/waffle_divorce_keyedarray.jl

@@ -0,0 +1,49 @@
+using StatsBase, NamedTupleTools
+using StanSample
+
+df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "WaffleDivorce.csv"), DataFrame);
+df.D = zscore(df.Divorce)
+df.M = zscore(df.Marriage)
+df.A = zscore(df.MedianAgeMarriage)
+data = (N=size(df, 1), D=df.D, A=df.A, M= df.M)
+
+stan5_1 = "
+data {
+    int < lower = 1 > N; // Sample size
+    vector[N] D; // Outcome
+    vector[N] A; // Predictor
+}
+parameters {
+    real a; // Intercept
+    real bA; // Slope (regression coefficients)
+    real < lower = 0 > sigma;    // Error SD
+}
+transformed parameters {
+    vector[N] mu;               // mu is a vector
+    for (i in 1:N)
+        mu[i] = a + bA * A[i];
+}
+model {
+    a ~ normal(0, 0.2);         //Priors
+    bA ~ normal(0, 0.5);
+    sigma ~ exponential(1);
+    D ~ normal(mu , sigma);     // Likelihood
+}
+generated quantities {
+    vector[N] log_lik;
+    for (i in 1:N)
+        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
+}
+";
+
+m5_1s = SampleModel("m5.1s", stan5_1)
+rc5_1s = stan_sample(m5_1s; data)
+
+if success(rc5_1s)
+
+    nt5_1s = read_samples(m5_1s, :particles)
+    NamedTupleTools.select(nt5_1s, (:a, :bA, :sigma)) |> display
+
+    chns = read_samples(m5_1s)
+    chns |> display
+end