From bc260b2d01274ceddb482a05056d537ae766f409 Mon Sep 17 00:00:00 2001 From: Rik Huijzer Date: Fri, 8 Jan 2021 22:43:28 +0100 Subject: [PATCH] Convert all existing models to Literate (#59) --- .github/workflows/CI.yml | 28 ++- Project.toml | 8 +- _css/jtd.css | 2 + index.md | 2 + models.md | 16 +- models/admit-reject.md | 35 +--- models/beta-binomial.md | 61 +----- models/chimpanzees.md | 37 +--- models/estimate-handedness-chimpanzees.md | 50 +---- models/ignoring-gender-admit.md | 1 - models/multinomial-poisson.md | 81 +------- models/multivariate-chimpanzees-priors.md | 207 ++------------------- models/non-centered-chimpanzees.md | 3 + models/non-identifiable.md | 42 +---- models/oceanic-tool-complexity.md | 46 +---- models/spatial-autocorrelation-oceanic.md | 70 +------ models/varying-intercepts-admission.md | 93 +-------- models/varying-intercepts-reedfrogs.md | 93 +-------- models/varying-slopes-cafe.md | 119 +----------- models/weakly-informative-priors.md | 35 +--- models/wild-chain.md | 24 +-- scripts/admit-reject.jl | 40 ++++ scripts/beta-binomial.jl | 46 +++++ scripts/chimpanzees.jl | 40 ++++ scripts/estimate-handedness-chimpanzees.jl | 57 ++++++ scripts/multinomial-poisson.jl | 70 +++++++ scripts/multivariate-chimpanzees-priors.jl | 207 +++++++++++++++++++++ scripts/non-identifiable.jl | 36 ++++ scripts/oceanic-tool-complexity.jl | 56 ++++++ scripts/spatial-autocorrelation-oceanic.jl | 73 ++++++++ scripts/varying-intercepts-admission.jl | 98 ++++++++++ scripts/varying-intercepts-reedfrogs.jl | 84 +++++++++ scripts/varying-slopes-cafe.jl | 109 +++++++++++ scripts/weakly-informative-priors.jl | 29 +++ scripts/wild-chain.jl | 21 +++ 35 files changed, 1056 insertions(+), 963 deletions(-) create mode 100644 scripts/admit-reject.jl create mode 100644 scripts/beta-binomial.jl create mode 100644 scripts/chimpanzees.jl create mode 100644 scripts/estimate-handedness-chimpanzees.jl create mode 100644 scripts/multinomial-poisson.jl create mode 100644 scripts/multivariate-chimpanzees-priors.jl create mode 100644 scripts/non-identifiable.jl create mode 100644 scripts/oceanic-tool-complexity.jl create mode 100644 scripts/spatial-autocorrelation-oceanic.jl create mode 100644 scripts/varying-intercepts-admission.jl create mode 100644 scripts/varying-intercepts-reedfrogs.jl create mode 100644 scripts/varying-slopes-cafe.jl create mode 100644 scripts/weakly-informative-priors.jl create mode 100644 scripts/wild-chain.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d5013d32..c83302fb 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -23,18 +23,30 @@ jobs: version: 1 # Latest stable release, for example, `1.5.3`. - name: Build site - run: julia -e ' - using Pkg; Pkg.add(["NodeJS", "Franklin"]); - using NodeJS; run(`$(npm_cmd()) install highlight.js`); - using Franklin; - Pkg.develop(path="."); Pkg.activate("."); Pkg.instantiate(); - Franklin.optimize(); - verify_links()' > build.log + run: | + julia -e ' + using Pkg; Pkg.add(["NodeJS", "Franklin"]); + using NodeJS; run(`$(npm_cmd()) install highlight.js`); + using Franklin; + Pkg.develop(path="."); Pkg.activate("."); Pkg.instantiate(); + Franklin.optimize(); + verify_links()' > build.log; + cat build.log - name: Deploy to secondary branch - # Only deploy when pushed to master. if: ${{ github.event_name == 'push' }} uses: peaceiris/actions-gh-pages@v3 with: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./__site + + - name: Validate output + run: | + # A Franklin Warning means an error occurred when evaluating code. + if grep -1 "Franklin Warning" build.log; then + echo "Franklin reported a warning" + exit 1 + else + echo "Franklin did not report a warning" + fi + diff --git a/Project.toml b/Project.toml index 4a77041d..5398f340 100644 --- a/Project.toml +++ b/Project.toml @@ -6,8 +6,13 @@ version = "1.1.3" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Franklin = "713c75ef-9fc9-4b05-94a9-213340da978e" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NodeJS = "2bd173c7-0d6d-553b-b6af-13a54713934c" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +RData = "df47a6cb-8c03-5eed-afd8-b6050d6c41da" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" @@ -15,9 +20,10 @@ Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" [compat] CSV = "0.8" DataFrames = "0.22" +Distributions = "0.24" Franklin = "0.10" NodeJS = "1.1" +RData = "0.7" StatsFuns = "0.9" StatsPlots = "0.14" Turing = "0.15" - diff --git a/_css/jtd.css b/_css/jtd.css index 2d2d04a4..1348b810 100644 --- a/_css/jtd.css +++ b/_css/jtd.css @@ -617,6 +617,7 @@ body { counter-reset: eqnum; } display: inline-block; } /* overwrite katex settings */ +/* Not interesting for this site. .katex-display::after { counter-increment: eqnum; content: "(" counter(eqnum) ")"; @@ -624,6 +625,7 @@ body { counter-reset: eqnum; } float: right; padding-right: 5px; } +*/ /* ================================================================== CODE & HIGHLIGHT.JS diff --git a/index.md b/index.md index 30b287e0..2450df45 100644 --- a/index.md +++ b/index.md @@ -11,6 +11,8 @@ We have tried our best to let this site be without mistakes. However, if you find a mistake on this website, please let us know. The code for this website is available on [GitHub](https://github.com/StatisticalRethinkingJulia/TuringModels.jl). +\toc + ### About the pages Each page aims to contain all the code required to reproduce the results. diff --git a/models.md b/models.md index 1ad240e3..d2b04793 100644 --- a/models.md +++ b/models.md @@ -21,7 +21,7 @@ Before you look at the models below, you might want to look at the [Basic Exampl - [m13.1](varying-intercepts-reedfrogs): Varying intercepts Reedfrogs - [m13.1](varying-slopes-cafe): Varying slopes cafes - [m13.2](multinomial-poisson): Multinomial Poisson regression -- [m13.3](varying-intercepts-admission): Varying intercepts admission decisions [not-implemented] +- [m13.3](varying-intercepts-admission): Varying intercepts admission decisions ## 1st Edition (2015) @@ -31,15 +31,15 @@ Before you look at the models below, you might want to look at the [Basic Exampl - [m8.2](wild-chain): Wild chain - [m8.3](weakly-informative-priors): Weakly informative priors - [m8.4](non-identifiable): Non-identifiable model -- [m10.3](chimpanzees): Chimpanzees [not-implemented] -- [m10.4](estimate-handedness-chimpanzees): Estimate handedness for each Chimpanzee [not-implemented] -- [m10.10](oceanic-tool-complexity): Oceanic tool complexity [not-implemented] -- [m10.yyt](admit-reject): Admit or reject [not-implemented] +- [m10.3](chimpanzees): Chimpanzees +- [m10.4](estimate-handedness-chimpanzees): Estimate handedness for each Chimpanzee +- [m10.10](oceanic-tool-complexity): Oceanic tool complexity +- [m10.yyt](admit-reject): Admit or reject - [m11.5](beta-binomial): Beta-binomial - [m11.7](multinomial-poisson): Multinomial Poisson regression - [m12.1](varying-intercepts-reedfrogs): Varying intercepts Reedfrogs - [m13.4](ignoring-gender-admit): Ignoring gender for admittance -- [m13.6](multivariate-chimpanzees-priors): Multivariate Chimpanzees priors [not-implemented] -- [m13.6nc](non-centered-chimpanzees): Non-centered Chimpanzees [not-implemented] -- [m13.7](spatial-autocorrelation-oceanic): Spatial autocorrelation in Oceanic tools [not-implemented] +- [m13.6](multivariate-chimpanzees-priors): Multivariate Chimpanzees priors +- [m13.6nc](non-centered-chimpanzees): Non-centered Chimpanzees +- [m13.7](spatial-autocorrelation-oceanic): Spatial autocorrelation in Oceanic tools - [m14.1](varying-slopes-cafe): Varying slopes cafes diff --git a/models/admit-reject.md b/models/admit-reject.md index d17a1269..5ea9ad5c 100644 --- a/models/admit-reject.md +++ b/models/admit-reject.md @@ -1,37 +1,6 @@ +++ title = "Admit or reject" +showall = true +++ -``` -using TuringModels, StatsFuns - -delim = ';' -d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "UCBadmit.csv"), DataFrame; delim); -size(d) # Should be 12x5 - -@model m_pois(admit, reject) = begin - α₁ ~ Normal(0,100) - α₂ ~ Normal(0,100) - - for i ∈ 1:length(admit) - λₐ = exp(α₁) - λᵣ = exp(α₂) - admit[i] ~ Poisson(λₐ) - reject[i] ~ Poisson(λᵣ) - end -end; - -chns = sample(m_pois(d[:, :admit], d[:, :reject]), Turing.NUTS(0.65), 1000); - -# Rethinking/CmdStan result - -m_10_yyt_result = " - mean sd 5.5% 94.5% n_eff Rhat - a1 4.99 0.02 4.95 5.02 2201 1 - a2 5.44 0.02 5.41 5.47 2468 1 -"; - -# Describe the draws - -chns |> display -``` +\literate{/scripts/admit-reject.jl} diff --git a/models/beta-binomial.md b/models/beta-binomial.md index 560a50a2..f2fb23e3 100644 --- a/models/beta-binomial.md +++ b/models/beta-binomial.md @@ -1,13 +1,15 @@ +++ title = "Beta-binomial" +showall = true +++ -In Rethinking 2nd Edition, this model is defined as +This beta-binomial model is `m11.5` in Statistical Rethinking 1st Edition. +In Statistical Rethinking 2nd Edition it is `m12.1` and defined as $$ \begin{aligned} A_i &\sim \text{BetaBinomial}(N_i, \overline{p}_i, \theta) \\ - \text{logit}(\overline{p}_i) &= \alpha_\text{GID}[i] \\ + \text{logit}(\overline{p}_i) &= \alpha_{\text{GID}[i]} \\ \alpha_j &\sim \text{Normal}(0, 1.5) \\ \theta &= \phi + 2 \\ \phi &\sim \text{Exponential(1)} @@ -16,58 +18,5 @@ $$ \toc -## Data +\literate{/scripts/beta-binomial.jl} -```julia:data -using DataFrames -using TuringModels # hide - -import CSV - -data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") -df = CSV.read(data_path, DataFrame; delim=';') - -@assert size(df) == (12, 5) # hide -write_csv(name, data) = CSV.write(joinpath(@OUTPUT, "$name.csv"), data) # hide -write_csv("data", df) # hide -``` -\output{data} - -\tableinput{}{./code/output/data.csv} - -## Model -```julia:model -using StatsFuns: logistic -using Turing - -@model function m11_5(admit, applications) - θ ~ truncated(Exponential(1), 0, Inf) - α ~ Normal(0, 2) - - # alpha and beta for the BetaBinomial must be provided. - # The two parameterizations are related by - # alpha = prob * theta, and beta = (1-prob) * theta. - # See https://github.com/rmcelreath/rethinking/blob/master/man/dbetabinom.Rd - - prob = logistic(α) - alpha = prob * θ - beta = (1 - prob) * θ - admit .~ BetaBinomial.(applications, alpha, beta) -end - -model = m11_5(df.admit, df.applications) -chains = sample(model, NUTS(0.65), 1000) -``` -\output{model} - -## Output - -\defaultoutput{} - -## Original output - -``` - mean sd 5.5% 94.5% n_eff Rhat -theta 2.74 0.96 1.43 4.37 3583 1 -a -0.37 0.31 -0.87 0.12 3210 1 -``` diff --git a/models/chimpanzees.md b/models/chimpanzees.md index ba2a79dc..2f71b629 100644 --- a/models/chimpanzees.md +++ b/models/chimpanzees.md @@ -1,39 +1,10 @@ +++ title = "Chimpanzees" +showall = true +++ -``` -using TuringModels, StatsFuns +This is the third Chimpanzees model (`m10.3`) in Statistical Rethinking Edition 1. -delim=';' -d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), DataFrame; delim); -size(d) # Should be 504x8 +\toc -# pulled_left, condition, prosoc_left -@model m10_3(y, x₁, x₂) = begin - α ~ Normal(0, 10) - βp ~ Normal(0, 10) - βpC ~ Normal(0, 10) - - logits = α .+ (βp .+ βpC * x₁) .* x₂ - y .~ BinomialLogit.(1, logits) -end; - -chns = sample(m10_3(d[:,:pulled_left], d[:,:condition], d[:,:prosoc_left]), - Turing.NUTS(0.65), 2000); - -# Rethinking result - -m_10_03t_result = " - Mean StdDev lower 0.89 upper 0.89 n_eff Rhat - a 0.05 0.13 -0.15 0.25 3284 1 - bp 0.62 0.22 0.28 0.98 3032 1 - bpC -0.11 0.26 -0.53 0.29 3184 1 -"; - -# Describe the draws - -chns |> display - -# End of m10.03t.jl -``` +\literate{/scripts/chimpanzees.jl} diff --git a/models/estimate-handedness-chimpanzees.md b/models/estimate-handedness-chimpanzees.md index 6cf638b3..fa5cd223 100644 --- a/models/estimate-handedness-chimpanzees.md +++ b/models/estimate-handedness-chimpanzees.md @@ -1,53 +1,11 @@ +++ title = "Estimate handedness for each Chimpanzee" +showall = true +++ -``` -using TuringModels, StatsFuns +This is model `m10.4` in Statistical Rethinking Edition 1. -delim=';' -d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), DataFrame; delim); -size(d) # Should be 504x8 +\toc -# pulled_left, actors, condition, prosoc_left -@model m10_4(y, actors, x₁, x₂) = begin - # Number of unique actors in the data set - N_actor = length(unique(actors)) +\literate{/scripts/estimate-handedness-chimpanzees.jl} - # Set an TArray for the priors/param - α ~ filldist(Normal(0, 10), N_actor) - βp ~ Normal(0, 10) - βpC ~ Normal(0, 10) - - logits = α[actors] .+ (βp .+ βpC * x₁) .* x₂ - y .~ BinomialLogit.(1, logits) -end; - -chns = sample(m10_4(d[:,:pulled_left], d[:,:actor],d[:,:condition], - d[:,:prosoc_left]), Turing.NUTS(0.65), 1000); - -# Rethinking/CmdStan results - -m_10_04s_result = " -Iterations = 1:1000 -Thinning interval = 1 -Chains = 1,2,3,4 -Samples per chain = 1000 - -Empirical Posterior Estimates: - Mean SD Naive SE MCSE ESS -a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000 -a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000 -a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000 -a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000 -a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000 -a.6 0.21599365 0.26307574 0.0041595927 0.0045153523 1000 -a.7 1.81090866 0.39318577 0.0062168129 0.0071483527 1000 -bp 0.83979926 0.26284676 0.0041559722 0.0059795826 1000 -bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000 -"; - -# Describe the draws - -chns |> display -``` diff --git a/models/ignoring-gender-admit.md b/models/ignoring-gender-admit.md index 2a832204..793aa945 100644 --- a/models/ignoring-gender-admit.md +++ b/models/ignoring-gender-admit.md @@ -1,7 +1,6 @@ +++ title = "Ignoring gender for admittance" showall = true -reeval = true +++ \literate{/scripts/ignoring-gender-admit.jl} diff --git a/models/multinomial-poisson.md b/models/multinomial-poisson.md index 5f7aa7fe..605139e4 100644 --- a/models/multinomial-poisson.md +++ b/models/multinomial-poisson.md @@ -1,5 +1,6 @@ +++ title = "Multinomial Poisson regression" +showall = true +++ On page 399 of McElreath (2015), `m13.2` is defined as @@ -17,82 +18,4 @@ $$ \toc -## Data - -```julia:data -import CSV - -using DataFrames -using TuringModels - -data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") -df = CSV.read(data_path, DataFrame; delim=';') -write_csv(name, data) = CSV.write(joinpath(@OUTPUT, "$name.csv"), data) # hide -write_csv("data", df) # hide -``` -\output{data} -\tableinput{}{./code/output/data.csv} - -## Model - -```julia:model -using Turing - -dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept))) -df.male = [g == "male" ? 1 : 0 for g in df.gender] -df.dept_id = [dept_map[de] for de in df.dept] - -@model m13_2(applications, dept_id, male, admit) = begin - sigma_dept ~ truncated(Cauchy(0, 2), 0, Inf) - bm ~ Normal(0, 1) - a ~ Normal(0, 10) - a_dept ~ filldist(Normal(a, sigma_dept), 6) - - logit_p = a_dept[dept_id] + bm*male - - admit .~ BinomialLogit.(applications, logit_p) -end - -chains = sample( - m13_2(df.applications, df.dept_id, df.male, df.admit), - Turing.NUTS(0.65), - 1000 -) -``` -\output{model} - -## Output - -\defaultoutput{} - -## Original output - -```text -Inference for Stan model: 359c2483e3bdbf74fd0484be27c2909b. - 3 chains, each with iter=4500; warmup=500; thin=1; - post-warmup draws per chain=4000, total post-warmup draws=12000. - - mean se_mean sd 2.5% 25% 50% 75% 97.5% - a_dept[1] 0.67 0.00 0.10 0.48 0.61 0.67 0.74 0.87 - a_dept[2] 0.63 0.00 0.12 0.40 0.55 0.63 0.71 0.85 - a_dept[3] -0.58 0.00 0.08 -0.73 -0.63 -0.58 -0.53 -0.44 - a_dept[4] -0.62 0.00 0.09 -0.79 -0.67 -0.62 -0.56 -0.45 - a_dept[5] -1.06 0.00 0.10 -1.26 -1.13 -1.06 -0.99 -0.86 - a_dept[6] -2.61 0.00 0.16 -2.92 -2.71 -2.60 -2.50 -2.31 - a -0.59 0.01 0.66 -1.94 -0.97 -0.59 -0.21 0.73 - bm -0.09 0.00 0.08 -0.25 -0.15 -0.09 -0.04 0.06 - sigma_dept 1.49 0.01 0.62 0.78 1.09 1.35 1.71 3.01 - lp__ -2602.27 0.04 2.27 -2607.74 -2603.48 -2601.90 -2600.61 -2598.95 - n_eff Rhat - a_dept[1] 8191 1 - a_dept[2] 8282 1 - a_dept[3] 10974 1 - a_dept[4] 10034 1 - a_dept[5] 10783 1 - a_dept[6] 10344 1 - a 5839 1 - bm 7137 1 - sigma_dept 4626 1 - lp__ 4200 1 -``` - +\literate{/scripts/multinomial-poisson.jl} diff --git a/models/multivariate-chimpanzees-priors.md b/models/multivariate-chimpanzees-priors.md index 0327bd4a..bc813962 100644 --- a/models/multivariate-chimpanzees-priors.md +++ b/models/multivariate-chimpanzees-priors.md @@ -1,206 +1,25 @@ +++ title = "Multivariate Chimpanzees priors" +showall = true +++ This is `m13.6` in Rethinking 1st Edition. In the 2nd Edition, the most similar one appears to be `m13.4nc`. +This model is not ran when generating these pages because the running time is about 20 minutes. +You can run it locally and inspect the output by downloading the file and running +``` +julia -i multivariate-chimpanzees.priors.jl ``` -using TuringModels - -# This script requires latest LKJ bijectors support. -# `] add Bijectors#master` to get latest Bijectors. - -data_path = joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv") -delim = ";" -d = CSV.read(data_path, DataFrame; delim) - -d.block_id = d.block - -@model m13_6(actor, block_id, condition, prosoc_left, pulled_left) = begin - # fixed priors - Rho_block ~ LKJ(3, 4.) - Rho_actor ~ LKJ(3, 4.) - sigma_block ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3) - sigma_actor ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3) - a ~ Normal(0, 1) - bp ~ Normal(0, 1) - bpc ~ Normal(0, 1) - - # adaptive priors - Sigma_block = sigma_block .* Rho_block .* sigma_block' - Sigma_block = (Sigma_block' + Sigma_block) / 2 - Sigma_actor = sigma_actor .* Rho_actor .* sigma_actor' - Sigma_actor = (Sigma_actor' + Sigma_actor) / 2 - - a_block_bp_block_bpc_block ~ filldist(MvNormal(zeros(3), Sigma_block), 6) - a_actor_bp_actor_bpc_actor ~ filldist(MvNormal(zeros(3), Sigma_actor), 7) - - a_block = a_block_bp_block_bpc_block[1, :] - bp_block = a_block_bp_block_bpc_block[2, :] - bpc_block = a_block_bp_block_bpc_block[3, :] - a_actor = a_actor_bp_actor_bpc_actor[1, :] - bp_actor = a_actor_bp_actor_bpc_actor[2, :] - bpc_actor = a_actor_bp_actor_bpc_actor[3, :] - - # linear models - BPC = bpc .+ bpc_actor[actor] + bpc_block[block_id] - BP = bp .+ bp_actor[actor] + bp_block[block_id] - A = a .+ a_actor[actor] + a_block[block_id] - logit_p = A + (BP + BPC .* condition) .* prosoc_left - - # likelihood - pulled_left .~ BinomialLogit.(1, logit_p) -end -chns = sample( - m13_6(d.actor, d.block_id, d.condition, d.prosoc_left, d.pulled_left), - Turing.NUTS(0.95), - 1000 -) +The file contains -chns |> display +```julia:code +# hideall +import TuringModels -m_13_6_rethinking = """ -Inference for Stan model: cdd1241666414818ec292db21291c409. - 3 chains, each with iter=5000; warmup=1000; thin=1; - post-warmup draws per chain=4000, total post-warmup draws=12000. - - mean se_mean sd 2.5% 25% 50% 75% 97.5% - bpc_actor[1] 0.02 0.01 0.41 -0.89 -0.16 0.01 0.20 0.91 - bpc_actor[2] 0.17 0.01 0.77 -1.20 -0.15 0.04 0.40 2.10 - bpc_actor[3] -0.28 0.01 0.48 -1.53 -0.49 -0.15 0.02 0.43 - bpc_actor[4] -0.05 0.01 0.43 -1.06 -0.22 -0.02 0.14 0.79 - bpc_actor[5] -0.07 0.01 0.43 -1.05 -0.25 -0.03 0.12 0.75 - bpc_actor[6] -0.12 0.01 0.44 -1.19 -0.30 -0.05 0.10 0.69 - bpc_actor[7] 0.40 0.02 0.70 -0.42 0.00 0.18 0.63 2.18 - bp_actor[1] 0.07 0.01 0.38 -0.67 -0.12 0.03 0.25 0.93 - bp_actor[2] 0.03 0.02 0.69 -1.22 -0.26 -0.01 0.26 1.54 - bp_actor[3] 0.24 0.01 0.42 -0.42 -0.02 0.14 0.44 1.29 - bp_actor[4] 0.17 0.01 0.40 -0.55 -0.05 0.09 0.36 1.11 - bp_actor[5] 0.10 0.01 0.38 -0.66 -0.09 0.05 0.28 0.99 - bp_actor[6] -0.32 0.01 0.43 -1.38 -0.56 -0.22 -0.02 0.30 - bp_actor[7] -0.04 0.01 0.45 -1.02 -0.24 -0.01 0.17 0.88 - a_actor[1] -1.04 0.02 0.69 -2.41 -1.50 -1.04 -0.58 0.30 - a_actor[2] 4.37 0.03 1.54 2.13 3.32 4.13 5.13 8.02 - a_actor[3] -1.33 0.02 0.70 -2.70 -1.79 -1.32 -0.85 0.00 - a_actor[4] -1.37 0.02 0.69 -2.75 -1.82 -1.35 -0.90 -0.04 - a_actor[5] -1.02 0.02 0.69 -2.38 -1.48 -1.01 -0.55 0.31 - a_actor[6] 0.18 0.02 0.69 -1.15 -0.29 0.18 0.62 1.54 - a_actor[7] 1.53 0.02 0.73 0.14 1.04 1.52 2.01 3.01 - bpc_block[1] -0.39 0.01 0.54 -1.69 -0.67 -0.25 -0.01 0.36 - bpc_block[2] 0.19 0.01 0.45 -0.57 -0.05 0.10 0.40 1.28 - bpc_block[3] 0.23 0.01 0.47 -0.52 -0.04 0.14 0.45 1.39 - bpc_block[4] -0.08 0.01 0.46 -1.14 -0.28 -0.04 0.13 0.83 - bpc_block[5] -0.08 0.01 0.42 -1.05 -0.28 -0.04 0.13 0.74 - bpc_block[6] 0.09 0.01 0.47 -0.86 -0.14 0.04 0.32 1.18 - bp_block[1] -0.18 0.01 0.44 -1.14 -0.42 -0.13 0.06 0.70 - bp_block[2] -0.12 0.01 0.42 -1.07 -0.33 -0.07 0.11 0.69 - bp_block[3] -0.24 0.01 0.45 -1.27 -0.48 -0.17 0.02 0.52 - bp_block[4] 0.16 0.01 0.43 -0.62 -0.08 0.10 0.37 1.15 - bp_block[5] 0.04 0.01 0.41 -0.78 -0.19 0.02 0.25 0.96 - bp_block[6] 0.64 0.02 0.57 -0.14 0.19 0.55 0.98 1.97 - a_block[1] -0.12 0.01 0.24 -0.71 -0.23 -0.07 0.01 0.24 - a_block[2] 0.06 0.00 0.22 -0.33 -0.04 0.03 0.16 0.61 - a_block[3] 0.10 0.01 0.24 -0.27 -0.02 0.05 0.20 0.69 - a_block[4] 0.00 0.00 0.20 -0.46 -0.09 0.00 0.09 0.42 - a_block[5] -0.02 0.00 0.21 -0.50 -0.11 -0.01 0.07 0.40 - a_block[6] 0.00 0.00 0.22 -0.48 -0.09 0.00 0.10 0.46 - a 0.26 0.02 0.63 -0.97 -0.16 0.26 0.69 1.50 - bp 0.71 0.01 0.41 -0.14 0.46 0.73 0.97 1.48 - bpc -0.03 0.01 0.43 -0.85 -0.31 -0.04 0.22 0.87 - sigma_actor[1] 2.31 0.01 0.88 1.18 1.72 2.14 2.69 4.52 - sigma_actor[2] 0.48 0.01 0.37 0.03 0.21 0.40 0.65 1.36 - sigma_actor[3] 0.51 0.01 0.45 0.03 0.19 0.40 0.71 1.69 - sigma_block[1] 0.24 0.01 0.21 0.01 0.09 0.19 0.33 0.77 - sigma_block[2] 0.59 0.01 0.41 0.04 0.30 0.51 0.79 1.58 - sigma_block[3] 0.53 0.01 0.41 0.03 0.23 0.44 0.73 1.55 - Rho_actor[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 - Rho_actor[1,2] -0.07 0.01 0.32 -0.66 -0.30 -0.07 0.16 0.54 - Rho_actor[1,3] 0.08 0.01 0.32 -0.56 -0.15 0.09 0.32 0.67 - Rho_actor[2,1] -0.07 0.01 0.32 -0.66 -0.30 -0.07 0.16 0.54 - Rho_actor[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 - Rho_actor[2,3] -0.06 0.00 0.31 -0.64 -0.29 -0.06 0.17 0.56 - Rho_actor[3,1] 0.08 0.01 0.32 -0.56 -0.15 0.09 0.32 0.67 - Rho_actor[3,2] -0.06 0.00 0.31 -0.64 -0.29 -0.06 0.17 0.56 - Rho_actor[3,3] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 - Rho_block[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 - Rho_block[1,2] -0.06 0.01 0.31 -0.65 -0.29 -0.07 0.16 0.56 - Rho_block[1,3] 0.03 0.00 0.32 -0.57 -0.20 0.03 0.26 0.63 - Rho_block[2,1] -0.06 0.01 0.31 -0.65 -0.29 -0.07 0.16 0.56 - Rho_block[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 - Rho_block[2,3] -0.05 0.00 0.32 -0.64 -0.27 -0.05 0.18 0.58 - Rho_block[3,1] 0.03 0.00 0.32 -0.57 -0.20 0.03 0.26 0.63 - Rho_block[3,2] -0.05 0.00 0.32 -0.64 -0.27 -0.05 0.18 0.58 - Rho_block[3,3] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 - lp__ -251.19 0.64 12.28 -273.25 -259.63 -252.09 -243.59 -225.28 - n_eff Rhat - bpc_actor[1] 4220 1.00 - bpc_actor[2] 2668 1.00 - bpc_actor[3] 1698 1.00 - bpc_actor[4] 2979 1.00 - bpc_actor[5] 2697 1.00 - bpc_actor[6] 2768 1.00 - bpc_actor[7] 1539 1.00 - bp_actor[1] 2359 1.00 - bp_actor[2] 1549 1.00 - bp_actor[3] 1926 1.00 - bp_actor[4] 2439 1.00 - bp_actor[5] 2525 1.00 - bp_actor[6] 2230 1.00 - bp_actor[7] 2745 1.00 - a_actor[1] 1355 1.00 - a_actor[2] 2638 1.00 - a_actor[3] 1439 1.00 - a_actor[4] 1413 1.00 - a_actor[5] 1337 1.00 - a_actor[6] 1345 1.00 - a_actor[7] 1554 1.00 - bpc_block[1] 1570 1.00 - bpc_block[2] 2224 1.00 - bpc_block[3] 2360 1.00 - bpc_block[4] 3057 1.00 - bpc_block[5] 2753 1.00 - bpc_block[6] 3163 1.00 - bp_block[1] 2282 1.00 - bp_block[2] 2251 1.00 - bp_block[3] 2085 1.00 - bp_block[4] 2411 1.00 - bp_block[5] 3092 1.00 - bp_block[6] 1413 1.00 - a_block[1] 1947 1.00 - a_block[2] 2493 1.00 - a_block[3] 1986 1.00 - a_block[4] 3390 1.00 - a_block[5] 3484 1.00 - a_block[6] 2141 1.00 - a 1313 1.00 - bp 2439 1.00 - bpc 2609 1.00 - sigma_actor[1] 3874 1.00 - sigma_actor[2] 814 1.00 - sigma_actor[3] 951 1.00 - sigma_block[1] 1096 1.00 - sigma_block[2] 1356 1.00 - sigma_block[3] 1109 1.00 - Rho_actor[1,1] NaN NaN - Rho_actor[1,2] 3766 1.00 - Rho_actor[1,3] 3052 1.00 - Rho_actor[2,1] 3766 1.00 - Rho_actor[2,2] 10943 1.00 - Rho_actor[2,3] 4210 1.00 - Rho_actor[3,1] 3052 1.00 - Rho_actor[3,2] 4210 1.00 - Rho_actor[3,3] 8994 1.00 - Rho_block[1,1] NaN NaN - Rho_block[1,2] 3441 1.00 - Rho_block[1,3] 5484 1.00 - Rho_block[2,1] 3441 1.00 - Rho_block[2,2] 11102 1.00 - Rho_block[2,3] 4267 1.00 - Rho_block[3,1] 5484 1.00 - Rho_block[3,2] 4267 1.00 - Rho_block[3,3] 7014 1.00 - lp__ 372 1.01 -""" +code_path = joinpath(TuringModels.project_root, "scripts", "multivariate-chimpanzees-priors.jl") +text = read(code_path, String) +print(text) ``` +\output{code} diff --git a/models/non-centered-chimpanzees.md b/models/non-centered-chimpanzees.md index 02e46de8..b0721d1d 100644 --- a/models/non-centered-chimpanzees.md +++ b/models/non-centered-chimpanzees.md @@ -2,6 +2,9 @@ title = "Non-centered Chimpanzees" +++ +This model is not executed on each change because it is similar to +[multivariate-chimpanzees-priors](/models/multivariate-chimpanzees-priors). + ``` using TuringModels using LinearAlgebra diff --git a/models/non-identifiable.md b/models/non-identifiable.md index f2d85283..ccc1b6a9 100644 --- a/models/non-identifiable.md +++ b/models/non-identifiable.md @@ -1,47 +1,9 @@ +++ title = "Non-identifiable model" +showall = true +++ \toc -## Data +\literate{/scripts/non-identifiable.jl} -```julia:data -using Distributions -using Random - -Random.seed!(1234) -y = rand(Normal(0,1), 100); -``` - -## Model - -```julia:model -using Turing - -# Can't really set a Uniform[-Inf,Inf] on σ - -@model m8_4(y) = begin - α₁ ~ Uniform(-3000, 1000) - α₂ ~ Uniform(-1000, 3000) - σ ~ truncated(Cauchy(0,1), 0, Inf) - - y .~ Normal(α₁ + α₂, σ) -end - -chains = sample(m8_4(y), NUTS(0.65), 2000) -``` -\output{model} - -## Output - -\defaultoutput{} - -## Original output - -``` - mean sd 5.5% 94.5% n_eff Rhat - a1 -861.15 558.17 -1841.89 -31.04 7 1.43 - a2 861.26 558.17 31.31 1842.00 7 1.43 - sigma 0.97 0.07 0.89 1.09 9 1.17 -``` diff --git a/models/oceanic-tool-complexity.md b/models/oceanic-tool-complexity.md index b11656cc..7daf891f 100644 --- a/models/oceanic-tool-complexity.md +++ b/models/oceanic-tool-complexity.md @@ -1,48 +1,8 @@ +++ title = "Oceanic tool complexity" +showall = true +++ -``` -using TuringModels +\toc -delim = ';' -d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "Kline.csv"), DataFrame; delim); -size(d) # Should be 10x5 - -d.log_pop = map(log, d.population) - -# New col contact_high, set binary values 1/0 if high/low contact -d.contact_high = [contact == "high" ? 1 : 0 for contact in d.contact] - -# This is supposed to be a "bad" model since we take non-centered data for the -# predictor log_pop -@model m10_10stan(total_tools, log_pop, contact_high) = begin - α ~ Normal(0, 100) - βp ~ Normal(0, 1) - βc ~ Normal(0, 1) - βpc ~ Normal(0, 1) - - for i ∈ 1:length(total_tools) - λ = exp(α + βp*log_pop[i] + βc*contact_high[i] + - βpc*contact_high[i]*log_pop[i]) - total_tools[i] ~ Poisson(λ) - end -end; - -posterior = sample(m10_10stan(d[:, :total_tools], d[:, :log_pop], - d[:, :contact_high]), Turing.NUTS(0.65), 1000); - -# Rethinking result - -m_10_10t_result = " - mean sd 5.5% 94.5% n_eff Rhat - a 0.94 0.37 0.36 1.53 3379 1 - bp 0.26 0.04 0.21 0.32 3347 1 - bc -0.08 0.84 -1.41 1.23 2950 1 - bpc 0.04 0.09 -0.10 0.19 2907 1 -"; - -# Describe the draws - -posterior |> display -``` +\literate{/scripts/oceanic-tool-complexity.jl} diff --git a/models/spatial-autocorrelation-oceanic.md b/models/spatial-autocorrelation-oceanic.md index 54b3863f..82f502d9 100644 --- a/models/spatial-autocorrelation-oceanic.md +++ b/models/spatial-autocorrelation-oceanic.md @@ -1,72 +1,8 @@ +++ title = "Spatial autocorrelation in Oceanic tools" +showall = true +++ -``` -using TuringModels -using RData -using LinearAlgebra +\toc -# This script requires latest LKJ bijectors support. -# `] add Bijectors#master` to get latest Bijectors. - -data_root = joinpath(@__DIR__, "..", "..", "data") -delim = ";" - -Kline2_path = joinpath(data_root, "Kline2.csv") -delim = ";" -d = CSV.read(Kline2_path, DataFrame; delim) - -Dmat_path = joinpath(data_root, "islandsDistMatrix.rda") -Dmat = load(Dmat_path)["islandsDistMatrix"] - -d.society = 1:10 - -@model m13_7(Dmat, society, logpop, total_tools) = begin - rhosq ~ truncated(Cauchy(0, 1), 0, Inf) - etasq ~ truncated(Cauchy(0, 1), 0, Inf) - bp ~ Normal(0, 1) - a ~ Normal(0, 10) - - # GPL2 - SIGMA_Dmat = etasq * exp.(-rhosq * Dmat.^2) - SIGMA_Dmat = SIGMA_Dmat + 0.01I - SIGMA_Dmat = (SIGMA_Dmat' + SIGMA_Dmat) / 2 - g ~ MvNormal(zeros(10), SIGMA_Dmat) - - log_lambda = a .+ g[society] .+ bp * logpop - - total_tools .~ Poisson.(exp.(log_lambda)) -end - -chns = sample( - m13_7(Dmat, d.society, d.logpop, d.total_tools), - Turing.NUTS(0.65), - 5000 -) - -chns |> display - -m_13_7_rethinking = """ -Inference for Stan model: 6422d8042e9cdd08dae2420ad26842f1. - 4 chains, each with iter=10000; warmup=2000; thin=1; - post-warmup draws per chain=8000, total post-warmup draws=32000. - - mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat - g[1] -0.27 0.01 0.44 -1.26 -0.50 -0.24 -0.01 0.53 3278 1 - g[2] -0.13 0.01 0.43 -1.07 -0.34 -0.10 0.12 0.67 3144 1 - g[3] -0.17 0.01 0.42 -1.10 -0.37 -0.14 0.07 0.59 3096 1 - g[4] 0.30 0.01 0.37 -0.51 0.12 0.30 0.50 1.02 3145 1 - g[5] 0.03 0.01 0.37 -0.77 -0.15 0.04 0.22 0.72 3055 1 - g[6] -0.46 0.01 0.38 -1.31 -0.64 -0.42 -0.23 0.19 3306 1 - g[7] 0.10 0.01 0.36 -0.70 -0.07 0.11 0.29 0.78 3175 1 - g[8] -0.26 0.01 0.37 -1.07 -0.43 -0.24 -0.06 0.40 3246 1 - g[9] 0.23 0.01 0.35 -0.52 0.07 0.25 0.42 0.89 3302 1 - g[10] -0.12 0.01 0.45 -1.04 -0.36 -0.11 0.13 0.77 5359 1 - a 1.31 0.02 1.15 -0.98 0.61 1.31 2.00 3.69 4389 1 - bp 0.25 0.00 0.11 0.02 0.18 0.24 0.31 0.47 5634 1 - etasq 0.34 0.01 0.53 0.03 0.10 0.20 0.39 1.52 5643 1 - rhosq 1.52 0.09 11.82 0.03 0.16 0.39 0.96 7.96 15955 1 - lp__ 925.98 0.03 2.96 919.16 924.20 926.34 928.14 930.67 7296 1 -""" -``` +\literate{/scripts/spatial-autocorrelation-oceanic.jl} diff --git a/models/varying-intercepts-admission.md b/models/varying-intercepts-admission.md index de2d334f..ab538002 100644 --- a/models/varying-intercepts-admission.md +++ b/models/varying-intercepts-admission.md @@ -1,95 +1,8 @@ +++ title = "Varying intercepts admission decisions" +showall = true +++ -``` -using TuringModels +\toc -# This script requires latest LKJ bijectors support. -# `] add Bijectors#master` to get latest Bijectors. - -data_path = joinpath(@__DIR__, "..", "..", "data", "UCBadmit.csv") -delim = ";" -d = CSV.read(data_path, DataFrame; delim) - -dept_map = Dict(key => idx for (idx, key) in enumerate(unique(d.dept))) -d.male = [g == "male" ? 1 : 0 for g in d.gender] -d.dept_id = [dept_map[de] for de in d.dept] - -@model m13_3(applications, dept_id, male, admit) = begin - Rho ~ LKJ(2, 2.) - sigma_dept ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2) - bm ~ Normal(0, 1) - a ~ Normal(0, 1) - - dist_mu = [a, bm] - dist_Sigma = sigma_dept .* Rho .* sigma_dept' - dist_Sigma = (dist_Sigma' + dist_Sigma) / 2 - - a_bm_dept ~ filldist(MvNormal(dist_mu, dist_Sigma), 6) - - a_dept, bm_dept = a_bm_dept[1, :], a_bm_dept[2, :] - logit_p = a_dept[dept_id] + bm_dept[dept_id] .* male - - admit .~ BinomialLogit.(applications, logit_p) -end - -chns = sample( - m13_3(d.applications, d.dept_id, d.male, d.admit), - Turing.NUTS(0.95), - 1000 -) - -m_13_3_rethinking = """ -Inference for Stan model: f0d86ec689cbf7921aab4fc0f55616d2. - 4 chains, each with iter=5000; warmup=1000; thin=1; - post-warmup draws per chain=4000, total post-warmup draws=16000. - - mean se_mean sd 2.5% 25% 50% 75% - bm_dept[1] -0.79 0.00 0.27 -1.32 -0.97 -0.79 -0.61 - bm_dept[2] -0.21 0.00 0.33 -0.87 -0.42 -0.20 0.00 - bm_dept[3] 0.08 0.00 0.14 -0.18 -0.01 0.08 0.18 - bm_dept[4] -0.09 0.00 0.14 -0.37 -0.19 -0.09 0.00 - bm_dept[5] 0.12 0.00 0.19 -0.24 -0.01 0.12 0.25 - bm_dept[6] -0.12 0.00 0.27 -0.67 -0.29 -0.12 0.05 - a_dept[1] 1.30 0.00 0.25 0.82 1.13 1.30 1.47 - a_dept[2] 0.74 0.00 0.33 0.08 0.52 0.73 0.95 - a_dept[3] -0.65 0.00 0.09 -0.82 -0.71 -0.65 -0.59 - a_dept[4] -0.62 0.00 0.10 -0.83 -0.69 -0.62 -0.55 - a_dept[5] -1.13 0.00 0.11 -1.36 -1.21 -1.13 -1.05 - a_dept[6] -2.60 0.00 0.20 -3.01 -2.73 -2.60 -2.46 - a -0.49 0.01 0.73 -1.96 -0.91 -0.49 -0.07 - bm -0.16 0.00 0.24 -0.65 -0.29 -0.16 -0.02 - sigma_dept[1] 1.67 0.01 0.63 0.88 1.25 1.53 1.94 - sigma_dept[2] 0.50 0.00 0.26 0.16 0.33 0.45 0.61 - Rho[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 - Rho[1,2] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08 - Rho[2,1] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08 - Rho[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 - lp__ -2593.92 0.04 3.17 -2601.03 -2595.85 -2593.57 -2591.65 - 97.5% n_eff Rhat - bm_dept[1] -0.27 6918 1 - bm_dept[2] 0.46 9133 1 - bm_dept[3] 0.36 13016 1 - bm_dept[4] 0.19 12258 1 - bm_dept[5] 0.49 12953 1 - bm_dept[6] 0.40 12011 1 - a_dept[1] 1.81 6993 1 - a_dept[2] 1.39 9335 1 - a_dept[3] -0.48 13489 1 - a_dept[4] -0.42 12595 1 - a_dept[5] -0.91 14659 1 - a_dept[6] -2.23 13168 1 - a 0.96 9541 1 - bm 0.31 8946 1 - sigma_dept[1] 3.27 9058 1 - sigma_dept[2] 1.15 6626 1 - Rho[1,1] 1.00 NaN NaN - Rho[1,2] 0.44 11358 1 - Rho[2,1] 0.44 11358 1 - Rho[2,2] 1.00 16069 1 - lp__ -2588.70 5310 1 -""" - -chns |> display -``` +\literate{/scripts/varying-intercepts-admission.jl} diff --git a/models/varying-intercepts-reedfrogs.md b/models/varying-intercepts-reedfrogs.md index d3e73f07..c86311ea 100644 --- a/models/varying-intercepts-reedfrogs.md +++ b/models/varying-intercepts-reedfrogs.md @@ -1,8 +1,10 @@ +++ title = "Varying intercepts Reedfrogs" +showall = true +++ On page 402 of Edition 2, this model is defined as + $$ \begin{aligned} S_i &\sim \text{Binomial}(N_i, p_i) \\ @@ -11,93 +13,6 @@ $$ \end{aligned} $$ -```julia:data -import CSV - -using DataFrames -using TuringModels - -data_path = joinpath(TuringModels.project_root, "data", "reedfrogs.csv") -df = CSV.read(data_path, DataFrame; delim=';'); -@assert size(df) == (48, 5) # hide -df.tank = 1:nrow(df) -write_csv(name, data) = CSV.write(joinpath(@OUTPUT, "$name.csv"), data) # hide -write_csv("data", df) # hide -``` -\output{data} - -\tableinput{}{./code/output/data.csv} - -```julia:model -using Turing - -@model reedfrogs(density, tank, surv, n) = begin - a_tank ~ filldist(Normal(0, 1.5), n) - - logitp = a_tank[tank] - surv .~ BinomialLogit.(density, logitp) -end - -n = nrow(df) -model = reedfrogs(df.density, df.tank, df.surv, n) -chains = sample(model, NUTS(0.65), 1000) -``` -\output{model} - -## Output - -\defaultoutput{} - -# Original output (Edition 1) +\toc -``` - mean sd 5.5% 94.5% n_eff Rhat -a_tank[1] 2.49 1.16 0.85 4.53 1079 1 -a_tank[2] 5.69 2.75 2.22 10.89 1055 1 -a_tank[3] 0.89 0.75 -0.23 2.16 1891 1 -a_tank[4] 5.71 2.70 2.21 10.85 684 1 -a_tank[5] 2.52 1.14 0.92 4.42 1640 1 -a_tank[6] 2.49 1.13 0.94 4.52 1164 1 -a_tank[7] 5.74 2.71 2.25 10.86 777 1 -a_tank[8] 2.52 1.19 0.95 4.42 1000 1 -a_tank[9] -0.46 0.69 -1.62 0.55 2673 1 -a_tank[10] 2.53 1.19 0.93 4.59 1430 1 -a_tank[11] 0.93 0.72 -0.17 2.11 1387 1 -a_tank[12] 0.47 0.74 -0.63 1.70 1346 1 -a_tank[13] 0.91 0.76 -0.25 2.30 1559 1 -a_tank[14] 0.00 0.66 -1.04 1.06 2085 1 -a_tank[15] 2.50 1.19 0.95 4.40 1317 1 -a_tank[16] 2.50 1.14 0.98 4.31 1412 1 -a_tank[17] 3.49 1.12 1.94 5.49 945 1 -a_tank[18] 2.59 0.75 1.50 3.81 1561 1 -a_tank[19] 2.11 0.64 1.15 3.15 1712 1 -a_tank[20] 6.40 2.57 3.11 11.04 996 1 -a_tank[21] 2.59 0.74 1.54 3.93 1233 1 -a_tank[22] 2.63 0.79 1.49 4.01 1184 1 -a_tank[23] 2.64 0.83 1.45 4.13 1379 1 -a_tank[24] 1.74 0.59 0.85 2.72 1736 1 -a_tank[25] -1.19 0.45 -1.90 -0.50 2145 1 -a_tank[26] 0.09 0.41 -0.53 0.78 2167 1 -a_tank[27] -1.75 0.56 -2.65 -0.88 1666 1 -a_tank[28] -0.58 0.43 -1.25 0.08 1567 1 -a_tank[29] 0.08 0.39 -0.54 0.71 3053 1 -a_tank[30] 1.43 0.49 0.66 2.24 2754 1 -a_tank[31] -0.79 0.44 -1.50 -0.12 1299 1 -a_tank[32] -0.42 0.41 -1.12 0.23 1661 1 -a_tank[33] 3.84 1.08 2.31 5.70 808 1 -a_tank[34] 3.00 0.85 1.83 4.36 1038 1 -a_tank[35] 2.96 0.82 1.82 4.25 1578 1 -a_tank[36] 2.14 0.55 1.31 3.08 1734 1 -a_tank[37] 2.12 0.56 1.31 3.04 1131 1 -a_tank[38] 6.72 2.62 3.45 11.44 706 1 -a_tank[39] 2.95 0.73 1.85 4.08 1509 1 -a_tank[40] 2.48 0.65 1.53 3.61 1731 1 -a_tank[41] -2.15 0.57 -3.11 -1.29 1231 1 -a_tank[42] -0.67 0.35 -1.22 -0.14 1444 1 -a_tank[43] -0.54 0.35 -1.12 0.03 1776 1 -a_tank[44] -0.43 0.34 -1.00 0.10 1735 1 -a_tank[45] 0.54 0.36 -0.04 1.14 1376 1 -a_tank[46] -0.67 0.34 -1.25 -0.15 1619 1 -a_tank[47] 2.14 0.55 1.31 3.04 1916 1 -a_tank[48] -0.06 0.35 -0.61 0.50 1932 1 -``` +\literate{/scripts/varying-intercepts-reedfrogs.jl} diff --git a/models/varying-slopes-cafe.md b/models/varying-slopes-cafe.md index 2400e13f..e1b1ba7d 100644 --- a/models/varying-slopes-cafe.md +++ b/models/varying-slopes-cafe.md @@ -19,121 +19,4 @@ $$ \toc -## Data - -```julia:data -import CSV - -using DataFrames -using TuringModels - -data_path = joinpath(TuringModels.project_root, "data", "d_13_1.csv") -df = CSV.read(data_path, DataFrame) -write_csv(name, data) = CSV.write(joinpath(@OUTPUT, "$name.csv"), data) # hide -write_csv("data", df) # hide -``` -\output{data} - -DataFrame `df` is shown in [Appendix I](#appendix_i). - -## Model - -```julia:model -using Turing - -@model m13_1(cafe, afternoon, wait) = begin - Rho ~ LKJ(2, 1.) - sigma ~ truncated(Cauchy(0, 2), 0, Inf) - sigma_cafe ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2) - a ~ Normal(0, 10) - b ~ Normal(0, 10) - - dist_mu = [a, b] - dist_Sigma = sigma_cafe .* Rho .* sigma_cafe' - dist_Sigma = (dist_Sigma' + dist_Sigma) / 2 - a_b_cafe ~ filldist(MvNormal(dist_mu, dist_Sigma), 20) - - a_cafe = a_b_cafe[1, :] - b_cafe = a_b_cafe[2, :] - - μ = a_cafe[cafe] + b_cafe[cafe] .* afternoon - wait .~ Normal.(μ, sigma) -end - -chains = sample( - m13_1(df.cafe, df.afternoon, df.wait), - # This model fails on NUTS(0.65). - Turing.NUTS(0.95), - 1000 -) -``` -\output{model} - -## Output - -\defaultoutput{} - -## Original output - -``` -Inference for Stan model: a73b0bd01032773825c6abf5575fd6e4. - 2 chains, each with iter=5000; warmup=2000; thin=1; - post-warmup draws per chain=3000, total post-warmup draws=6000. - - mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat - b_cafe[1] -1.29 0.00 0.18 -1.69 -1.40 -1.28 -1.18 -0.96 2548 1.00 - b_cafe[2] -1.20 0.00 0.18 -1.57 -1.31 -1.21 -1.09 -0.83 3288 1.00 - b_cafe[3] -1.26 0.00 0.18 -1.63 -1.36 -1.25 -1.15 -0.91 4302 1.00 - b_cafe[4] -1.29 0.00 0.18 -1.68 -1.39 -1.28 -1.18 -0.96 2960 1.00 - b_cafe[5] -1.26 0.00 0.20 -1.68 -1.38 -1.25 -1.14 -0.88 3406 1.00 - b_cafe[6] -1.28 0.00 0.18 -1.66 -1.38 -1.27 -1.17 -0.93 2983 1.00 - b_cafe[7] -1.23 0.00 0.18 -1.62 -1.33 -1.23 -1.13 -0.86 4430 1.00 - b_cafe[8] -1.25 0.00 0.17 -1.62 -1.35 -1.24 -1.15 -0.90 3838 1.00 - b_cafe[9] -1.13 0.01 0.19 -1.45 -1.26 -1.16 -1.02 -0.70 1372 1.00 - b_cafe[10] -1.19 0.00 0.18 -1.54 -1.30 -1.20 -1.09 -0.80 3924 1.00 - b_cafe[11] -1.04 0.01 0.22 -1.38 -1.20 -1.07 -0.90 -0.55 416 1.01 - b_cafe[12] -1.21 0.00 0.18 -1.55 -1.32 -1.22 -1.11 -0.82 3781 1.00 - b_cafe[13] -1.32 0.00 0.19 -1.76 -1.43 -1.30 -1.20 -0.99 1509 1.00 - b_cafe[14] -1.37 0.01 0.20 -1.82 -1.49 -1.34 -1.23 -1.04 760 1.01 - b_cafe[15] -1.52 0.02 0.27 -2.11 -1.70 -1.49 -1.31 -1.13 161 1.02 - b_cafe[16] -1.18 0.00 0.18 -1.51 -1.29 -1.20 -1.08 -0.79 3430 1.00 - b_cafe[17] -1.16 0.00 0.19 -1.50 -1.28 -1.18 -1.05 -0.73 2034 1.00 - b_cafe[18] -1.28 0.00 0.21 -1.72 -1.41 -1.28 -1.16 -0.86 3166 1.00 - b_cafe[19] -1.05 0.01 0.22 -1.38 -1.21 -1.09 -0.91 -0.55 363 1.01 - b_cafe[20] -1.10 0.01 0.20 -1.43 -1.23 -1.12 -0.98 -0.64 613 1.01 - a_cafe[1] 4.08 0.00 0.18 3.74 3.97 4.08 4.20 4.44 5407 1.00 - a_cafe[2] 2.38 0.00 0.18 2.03 2.26 2.38 2.50 2.72 5091 1.00 - a_cafe[3] 3.94 0.00 0.18 3.60 3.82 3.94 4.06 4.30 6412 1.00 - a_cafe[4] 3.45 0.00 0.18 3.11 3.34 3.45 3.57 3.81 5699 1.00 - a_cafe[5] 2.14 0.00 0.18 1.79 2.02 2.14 2.26 2.50 5380 1.00 - a_cafe[6] 4.26 0.00 0.17 3.92 4.15 4.26 4.38 4.61 5192 1.00 - a_cafe[7] 3.56 0.00 0.18 3.21 3.44 3.56 3.68 3.91 5495 1.00 - a_cafe[8] 3.79 0.00 0.18 3.44 3.68 3.79 3.91 4.14 5661 1.00 - a_cafe[9] 3.89 0.00 0.18 3.53 3.77 3.89 4.01 4.23 4135 1.00 - a_cafe[10] 3.69 0.00 0.18 3.34 3.57 3.69 3.81 4.05 5761 1.00 - a_cafe[11] 2.47 0.00 0.19 2.11 2.35 2.48 2.60 2.84 1548 1.00 - a_cafe[12] 4.08 0.00 0.17 3.74 3.97 4.08 4.20 4.42 5712 1.00 - a_cafe[13] 3.88 0.00 0.18 3.52 3.75 3.87 4.00 4.24 4061 1.00 - a_cafe[14] 3.33 0.00 0.19 2.97 3.21 3.33 3.46 3.70 3082 1.00 - a_cafe[15] 4.23 0.01 0.20 3.85 4.09 4.22 4.36 4.65 471 1.01 - a_cafe[16] 3.60 0.00 0.17 3.25 3.48 3.60 3.71 3.93 5527 1.00 - a_cafe[17] 4.43 0.00 0.18 4.09 4.31 4.43 4.55 4.78 5089 1.00 - a_cafe[18] 6.10 0.00 0.19 5.73 5.97 6.09 6.22 6.46 4780 1.00 - a_cafe[19] 3.50 0.00 0.19 3.12 3.38 3.50 3.63 3.86 1800 1.00 - a_cafe[20] 3.90 0.00 0.18 3.53 3.78 3.90 4.03 4.25 3465 1.00 - a 3.73 0.00 0.21 3.32 3.60 3.73 3.87 4.16 7131 1.00 - b -1.23 0.00 0.09 -1.40 -1.29 -1.23 -1.18 -1.06 2021 1.00 - sigma_cafe[1] 0.91 0.00 0.17 0.65 0.80 0.89 1.01 1.29 5579 1.00 - sigma_cafe[2] 0.21 0.01 0.12 0.01 0.12 0.20 0.29 0.46 72 1.05 - sigma 0.49 0.00 0.03 0.44 0.47 0.49 0.51 0.55 2271 1.00 - Rho[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN - Rho[1,2] -0.17 0.01 0.34 -0.75 -0.42 -0.20 0.06 0.57 3300 1.00 - Rho[2,1] -0.17 0.01 0.34 -0.75 -0.42 -0.20 0.06 0.57 3300 1.00 - Rho[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 5777 1.00 - lp__ 62.41 3.30 17.59 41.97 52.03 58.04 66.26 118.46 28 1.12 -``` - -## Appendix - -### Appendix I -\tableinput{}{./code/output/data.csv} +\literate{/scripts/varying-slopes-cafe.jl} diff --git a/models/weakly-informative-priors.md b/models/weakly-informative-priors.md index 32426c86..c4d771af 100644 --- a/models/weakly-informative-priors.md +++ b/models/weakly-informative-priors.md @@ -6,37 +6,4 @@ This model is a continuation of the [Wild Chain](/models/wild-chain). \toc -## Data - -```julia:data -y = [-1,1] -``` - -## Model - -```julia:model -using Turing - -@model function m8_3(y) - α ~ Normal(1, 10) - σ ~ truncated(Cauchy(0, 1), 0, Inf) - - μ = α - y .~ Normal(μ, σ) -end - -chains = sample(m8_3(y), NUTS(0.65), 1000); -``` -\output{model} - -## Output - -\defaultoutput{} - -## Original output - -``` - mean sd 5.5% 94.5% n_eff Rhat -alpha 0.09 1.63 -2.13 2.39 959 1 -sigma 2.04 2.05 0.68 4.83 1090 1 -``` +\literate{/scripts/weakly-informative-priors.jl} diff --git a/models/wild-chain.md b/models/wild-chain.md index ba821de8..1b784b03 100644 --- a/models/wild-chain.md +++ b/models/wild-chain.md @@ -6,26 +6,4 @@ This model shows what happens if you use extremely flat priors, and is fixed in \toc -## Data - -```julia:data -y = [-1, 1] -``` - -## Model -```julia:model -using Turing - -@model m8_2(y) = begin - σ ~ FlatPos(0.0) # improper prior with probability one everywhere above 0.0 - α ~ Flat() # improper prior with pobability one everywhere - - y .~ Normal(α, σ) -end - -chains = sample(m8_2(y), NUTS(0.65), 1000) -``` - -## Output - -\defaultoutput{} +\literate{/scripts/wild-chain.jl} diff --git a/scripts/admit-reject.jl b/scripts/admit-reject.jl new file mode 100644 index 00000000..1ad568e5 --- /dev/null +++ b/scripts/admit-reject.jl @@ -0,0 +1,40 @@ +# ## Data + +using CSV +using DataFrames +using StatsFuns +using Turing + +import TuringModels + +file_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") +df = CSV.read(file_path, DataFrame; delim=';') + +# ## Model + +@model m_pois(admit, reject) = begin + α₁ ~ Normal(0,100) + α₂ ~ Normal(0,100) + + for i ∈ 1:length(admit) + λₐ = exp(α₁) + λᵣ = exp(α₂) + admit[i] ~ Poisson(λₐ) + reject[i] ~ Poisson(λᵣ) + end +end; + +# ## Output + +chains = sample(m_pois(df.admit, df.reject), NUTS(0.65), 1000) + +# \defaultoutput{} + +# ## Original output + +m_10_yyt_result = " + mean sd 5.5% 94.5% n_eff Rhat + a1 4.99 0.02 4.95 5.02 2201 1 + a2 5.44 0.02 5.41 5.47 2468 1 +"; + diff --git a/scripts/beta-binomial.jl b/scripts/beta-binomial.jl new file mode 100644 index 00000000..55678780 --- /dev/null +++ b/scripts/beta-binomial.jl @@ -0,0 +1,46 @@ + +# ## Data + +using DataFrames + +import CSV +import TuringModels + +data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") +df = CSV.read(data_path, DataFrame; delim=';') + +# ## Model + +using StatsFuns: logistic +using Turing + +@model function m11_5(admit, applications) + θ ~ truncated(Exponential(1), 0, Inf) + α ~ Normal(0, 2) + + ## alpha and beta for the BetaBinomial must be provided. + ## The two parameterizations are related by + ## alpha = prob * theta, and beta = (1-prob) * theta. + ## See https://github.com/rmcelreath/rethinking/blob/master/man/dbetabinom.Rd + + prob = logistic(α) + alpha = prob * θ + beta = (1 - prob) * θ + admit .~ BetaBinomial.(applications, alpha, beta) +end; + +model = m11_5(df.admit, df.applications) + +# ## Output + +chains = sample(model, NUTS(0.65), 1000) + +# \defaultoutput{} + +# ## Original output + +""" + mean sd 5.5% 94.5% n_eff Rhat +theta 2.74 0.96 1.43 4.37 3583 1 +a -0.37 0.31 -0.87 0.12 3210 1 +""" diff --git a/scripts/chimpanzees.jl b/scripts/chimpanzees.jl new file mode 100644 index 00000000..3a1a6101 --- /dev/null +++ b/scripts/chimpanzees.jl @@ -0,0 +1,40 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames + +file_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv") +df = CSV.read(file_path, DataFrame; delim=';'); + +# ## Model + +using StatsFuns +using Turing + +@model function m10_3(y, x₁, x₂) + α ~ Normal(0, 10) + βp ~ Normal(0, 10) + βpC ~ Normal(0, 10) + + logits = α .+ (βp .+ βpC * x₁) .* x₂ + y .~ BinomialLogit.(1, logits) +end +model = m10_3(df.pulled_left, df.condition, df.prosoc_left) + +# ## Output + +chains = sample(model, NUTS(0.65), 2000) + +# \defaultoutput{} + +# ## Original output + +""" + Mean StdDev lower 0.89 upper 0.89 n_eff Rhat + a 0.05 0.13 -0.15 0.25 3284 1 + bp 0.62 0.22 0.28 0.98 3032 1 + bpC -0.11 0.26 -0.53 0.29 3184 1 +"""; + diff --git a/scripts/estimate-handedness-chimpanzees.jl b/scripts/estimate-handedness-chimpanzees.jl new file mode 100644 index 00000000..24126a49 --- /dev/null +++ b/scripts/estimate-handedness-chimpanzees.jl @@ -0,0 +1,57 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames +using StatsFuns + +data_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv") +df = CSV.read(data_path, DataFrame; delim=';') + +# ## Model + +using Turing + +@model m10_4(y, actors, x₁, x₂) = begin + ## Number of unique actors in the data set + N_actor = length(unique(actors)) + + ## Set an TArray for the priors/param + α ~ filldist(Normal(0, 10), N_actor) + βp ~ Normal(0, 10) + βpC ~ Normal(0, 10) + + logits = α[actors] .+ (βp .+ βpC * x₁) .* x₂ + y .~ BinomialLogit.(1, logits) +end + +model = m10_4(df.pulled_left, df.actor, df.condition, df.prosoc_left); + +# ## Output + +chains = sample(model, NUTS(0.65), 1000) + +# \defaultoutput{} + +# ## Original output + +m_10_04s_result = " +Iterations = 1:1000 +Thinning interval = 1 +Chains = 1,2,3,4 +Samples per chain = 1000 + +Empirical Posterior Estimates: + Mean SD Naive SE MCSE ESS +a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000 +a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000 +a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000 +a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000 +a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000 +a.6 0.21599365 0.26307574 0.0041595927 0.0045153523 1000 +a.7 1.81090866 0.39318577 0.0062168129 0.0071483527 1000 +bp 0.83979926 0.26284676 0.0041559722 0.0059795826 1000 +bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000 +"; + diff --git a/scripts/multinomial-poisson.jl b/scripts/multinomial-poisson.jl new file mode 100644 index 00000000..274b0f61 --- /dev/null +++ b/scripts/multinomial-poisson.jl @@ -0,0 +1,70 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames + +data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") +df = CSV.read(data_path, DataFrame; delim=';') + +# ## Model + +using Turing + +dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept))) +df.male = [g == "male" ? 1 : 0 for g in df.gender] +df.dept_id = [dept_map[de] for de in df.dept] + +@model m13_2(applications, dept_id, male, admit) = begin + sigma_dept ~ truncated(Cauchy(0, 2), 0, Inf) + bm ~ Normal(0, 1) + a ~ Normal(0, 10) + a_dept ~ filldist(Normal(a, sigma_dept), 6) + + logit_p = a_dept[dept_id] + bm*male + + admit .~ BinomialLogit.(applications, logit_p) +end + +chains = sample( + m13_2(df.applications, df.dept_id, df.male, df.admit), + Turing.NUTS(0.65), + 1000 +) + +# ## Output + +# \defaultoutput{} + +# ## Original output + +""" +Inference for Stan model: 359c2483e3bdbf74fd0484be27c2909b. + 3 chains, each with iter=4500; warmup=500; thin=1; + post-warmup draws per chain=4000, total post-warmup draws=12000. + + mean se_mean sd 2.5% 25% 50% 75% 97.5% + a_dept[1] 0.67 0.00 0.10 0.48 0.61 0.67 0.74 0.87 + a_dept[2] 0.63 0.00 0.12 0.40 0.55 0.63 0.71 0.85 + a_dept[3] -0.58 0.00 0.08 -0.73 -0.63 -0.58 -0.53 -0.44 + a_dept[4] -0.62 0.00 0.09 -0.79 -0.67 -0.62 -0.56 -0.45 + a_dept[5] -1.06 0.00 0.10 -1.26 -1.13 -1.06 -0.99 -0.86 + a_dept[6] -2.61 0.00 0.16 -2.92 -2.71 -2.60 -2.50 -2.31 + a -0.59 0.01 0.66 -1.94 -0.97 -0.59 -0.21 0.73 + bm -0.09 0.00 0.08 -0.25 -0.15 -0.09 -0.04 0.06 + sigma_dept 1.49 0.01 0.62 0.78 1.09 1.35 1.71 3.01 + lp__ -2602.27 0.04 2.27 -2607.74 -2603.48 -2601.90 -2600.61 -2598.95 + n_eff Rhat + a_dept[1] 8191 1 + a_dept[2] 8282 1 + a_dept[3] 10974 1 + a_dept[4] 10034 1 + a_dept[5] 10783 1 + a_dept[6] 10344 1 + a 5839 1 + bm 7137 1 + sigma_dept 4626 1 + lp__ 4200 1 +"""; + diff --git a/scripts/multivariate-chimpanzees-priors.jl b/scripts/multivariate-chimpanzees-priors.jl new file mode 100644 index 00000000..2caacc29 --- /dev/null +++ b/scripts/multivariate-chimpanzees-priors.jl @@ -0,0 +1,207 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames + +data_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv") +df = CSV.read(data_path, DataFrame; delim=';') + +df.block_id = df.block + +# ## Model + +using Turing + +@model m13_6(actor, block_id, condition, prosoc_left, pulled_left) = begin + ## fixed priors + Rho_block ~ LKJ(3, 4.) + Rho_actor ~ LKJ(3, 4.) + sigma_block ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3) + sigma_actor ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3) + a ~ Normal(0, 1) + bp ~ Normal(0, 1) + bpc ~ Normal(0, 1) + + ## adaptive priors + Sigma_block = sigma_block .* Rho_block .* sigma_block' + Sigma_block = (Sigma_block' + Sigma_block) / 2 + Sigma_actor = sigma_actor .* Rho_actor .* sigma_actor' + Sigma_actor = (Sigma_actor' + Sigma_actor) / 2 + + a_block_bp_block_bpc_block ~ filldist(MvNormal(zeros(3), Sigma_block), 6) + a_actor_bp_actor_bpc_actor ~ filldist(MvNormal(zeros(3), Sigma_actor), 7) + + a_block = a_block_bp_block_bpc_block[1, :] + bp_block = a_block_bp_block_bpc_block[2, :] + bpc_block = a_block_bp_block_bpc_block[3, :] + a_actor = a_actor_bp_actor_bpc_actor[1, :] + bp_actor = a_actor_bp_actor_bpc_actor[2, :] + bpc_actor = a_actor_bp_actor_bpc_actor[3, :] + + ## linear models + BPC = bpc .+ bpc_actor[actor] + bpc_block[block_id] + BP = bp .+ bp_actor[actor] + bp_block[block_id] + A = a .+ a_actor[actor] + a_block[block_id] + logit_p = A + (BP + BPC .* condition) .* prosoc_left + + ## likelihood + pulled_left .~ BinomialLogit.(1, logit_p) +end + +# ## Output + +chains = sample( + m13_6(df.actor, df.block_id, df.condition, df.prosoc_left, df.pulled_left), + NUTS(0.95), + 1000 +) + +# \defaultoutput{} + +# ## Original output + +m_13_6_rethinking = """ +Inference for Stan model: cdd1241666414818ec292db21291c409. + 3 chains, each with iter=5000; warmup=1000; thin=1; + post-warmup draws per chain=4000, total post-warmup draws=12000. + + mean se_mean sd 2.5% 25% 50% 75% 97.5% + bpc_actor[1] 0.02 0.01 0.41 -0.89 -0.16 0.01 0.20 0.91 + bpc_actor[2] 0.17 0.01 0.77 -1.20 -0.15 0.04 0.40 2.10 + bpc_actor[3] -0.28 0.01 0.48 -1.53 -0.49 -0.15 0.02 0.43 + bpc_actor[4] -0.05 0.01 0.43 -1.06 -0.22 -0.02 0.14 0.79 + bpc_actor[5] -0.07 0.01 0.43 -1.05 -0.25 -0.03 0.12 0.75 + bpc_actor[6] -0.12 0.01 0.44 -1.19 -0.30 -0.05 0.10 0.69 + bpc_actor[7] 0.40 0.02 0.70 -0.42 0.00 0.18 0.63 2.18 + bp_actor[1] 0.07 0.01 0.38 -0.67 -0.12 0.03 0.25 0.93 + bp_actor[2] 0.03 0.02 0.69 -1.22 -0.26 -0.01 0.26 1.54 + bp_actor[3] 0.24 0.01 0.42 -0.42 -0.02 0.14 0.44 1.29 + bp_actor[4] 0.17 0.01 0.40 -0.55 -0.05 0.09 0.36 1.11 + bp_actor[5] 0.10 0.01 0.38 -0.66 -0.09 0.05 0.28 0.99 + bp_actor[6] -0.32 0.01 0.43 -1.38 -0.56 -0.22 -0.02 0.30 + bp_actor[7] -0.04 0.01 0.45 -1.02 -0.24 -0.01 0.17 0.88 + a_actor[1] -1.04 0.02 0.69 -2.41 -1.50 -1.04 -0.58 0.30 + a_actor[2] 4.37 0.03 1.54 2.13 3.32 4.13 5.13 8.02 + a_actor[3] -1.33 0.02 0.70 -2.70 -1.79 -1.32 -0.85 0.00 + a_actor[4] -1.37 0.02 0.69 -2.75 -1.82 -1.35 -0.90 -0.04 + a_actor[5] -1.02 0.02 0.69 -2.38 -1.48 -1.01 -0.55 0.31 + a_actor[6] 0.18 0.02 0.69 -1.15 -0.29 0.18 0.62 1.54 + a_actor[7] 1.53 0.02 0.73 0.14 1.04 1.52 2.01 3.01 + bpc_block[1] -0.39 0.01 0.54 -1.69 -0.67 -0.25 -0.01 0.36 + bpc_block[2] 0.19 0.01 0.45 -0.57 -0.05 0.10 0.40 1.28 + bpc_block[3] 0.23 0.01 0.47 -0.52 -0.04 0.14 0.45 1.39 + bpc_block[4] -0.08 0.01 0.46 -1.14 -0.28 -0.04 0.13 0.83 + bpc_block[5] -0.08 0.01 0.42 -1.05 -0.28 -0.04 0.13 0.74 + bpc_block[6] 0.09 0.01 0.47 -0.86 -0.14 0.04 0.32 1.18 + bp_block[1] -0.18 0.01 0.44 -1.14 -0.42 -0.13 0.06 0.70 + bp_block[2] -0.12 0.01 0.42 -1.07 -0.33 -0.07 0.11 0.69 + bp_block[3] -0.24 0.01 0.45 -1.27 -0.48 -0.17 0.02 0.52 + bp_block[4] 0.16 0.01 0.43 -0.62 -0.08 0.10 0.37 1.15 + bp_block[5] 0.04 0.01 0.41 -0.78 -0.19 0.02 0.25 0.96 + bp_block[6] 0.64 0.02 0.57 -0.14 0.19 0.55 0.98 1.97 + a_block[1] -0.12 0.01 0.24 -0.71 -0.23 -0.07 0.01 0.24 + a_block[2] 0.06 0.00 0.22 -0.33 -0.04 0.03 0.16 0.61 + a_block[3] 0.10 0.01 0.24 -0.27 -0.02 0.05 0.20 0.69 + a_block[4] 0.00 0.00 0.20 -0.46 -0.09 0.00 0.09 0.42 + a_block[5] -0.02 0.00 0.21 -0.50 -0.11 -0.01 0.07 0.40 + a_block[6] 0.00 0.00 0.22 -0.48 -0.09 0.00 0.10 0.46 + a 0.26 0.02 0.63 -0.97 -0.16 0.26 0.69 1.50 + bp 0.71 0.01 0.41 -0.14 0.46 0.73 0.97 1.48 + bpc -0.03 0.01 0.43 -0.85 -0.31 -0.04 0.22 0.87 + sigma_actor[1] 2.31 0.01 0.88 1.18 1.72 2.14 2.69 4.52 + sigma_actor[2] 0.48 0.01 0.37 0.03 0.21 0.40 0.65 1.36 + sigma_actor[3] 0.51 0.01 0.45 0.03 0.19 0.40 0.71 1.69 + sigma_block[1] 0.24 0.01 0.21 0.01 0.09 0.19 0.33 0.77 + sigma_block[2] 0.59 0.01 0.41 0.04 0.30 0.51 0.79 1.58 + sigma_block[3] 0.53 0.01 0.41 0.03 0.23 0.44 0.73 1.55 + Rho_actor[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 + Rho_actor[1,2] -0.07 0.01 0.32 -0.66 -0.30 -0.07 0.16 0.54 + Rho_actor[1,3] 0.08 0.01 0.32 -0.56 -0.15 0.09 0.32 0.67 + Rho_actor[2,1] -0.07 0.01 0.32 -0.66 -0.30 -0.07 0.16 0.54 + Rho_actor[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 + Rho_actor[2,3] -0.06 0.00 0.31 -0.64 -0.29 -0.06 0.17 0.56 + Rho_actor[3,1] 0.08 0.01 0.32 -0.56 -0.15 0.09 0.32 0.67 + Rho_actor[3,2] -0.06 0.00 0.31 -0.64 -0.29 -0.06 0.17 0.56 + Rho_actor[3,3] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 + Rho_block[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 + Rho_block[1,2] -0.06 0.01 0.31 -0.65 -0.29 -0.07 0.16 0.56 + Rho_block[1,3] 0.03 0.00 0.32 -0.57 -0.20 0.03 0.26 0.63 + Rho_block[2,1] -0.06 0.01 0.31 -0.65 -0.29 -0.07 0.16 0.56 + Rho_block[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 + Rho_block[2,3] -0.05 0.00 0.32 -0.64 -0.27 -0.05 0.18 0.58 + Rho_block[3,1] 0.03 0.00 0.32 -0.57 -0.20 0.03 0.26 0.63 + Rho_block[3,2] -0.05 0.00 0.32 -0.64 -0.27 -0.05 0.18 0.58 + Rho_block[3,3] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 + lp__ -251.19 0.64 12.28 -273.25 -259.63 -252.09 -243.59 -225.28 + n_eff Rhat + bpc_actor[1] 4220 1.00 + bpc_actor[2] 2668 1.00 + bpc_actor[3] 1698 1.00 + bpc_actor[4] 2979 1.00 + bpc_actor[5] 2697 1.00 + bpc_actor[6] 2768 1.00 + bpc_actor[7] 1539 1.00 + bp_actor[1] 2359 1.00 + bp_actor[2] 1549 1.00 + bp_actor[3] 1926 1.00 + bp_actor[4] 2439 1.00 + bp_actor[5] 2525 1.00 + bp_actor[6] 2230 1.00 + bp_actor[7] 2745 1.00 + a_actor[1] 1355 1.00 + a_actor[2] 2638 1.00 + a_actor[3] 1439 1.00 + a_actor[4] 1413 1.00 + a_actor[5] 1337 1.00 + a_actor[6] 1345 1.00 + a_actor[7] 1554 1.00 + bpc_block[1] 1570 1.00 + bpc_block[2] 2224 1.00 + bpc_block[3] 2360 1.00 + bpc_block[4] 3057 1.00 + bpc_block[5] 2753 1.00 + bpc_block[6] 3163 1.00 + bp_block[1] 2282 1.00 + bp_block[2] 2251 1.00 + bp_block[3] 2085 1.00 + bp_block[4] 2411 1.00 + bp_block[5] 3092 1.00 + bp_block[6] 1413 1.00 + a_block[1] 1947 1.00 + a_block[2] 2493 1.00 + a_block[3] 1986 1.00 + a_block[4] 3390 1.00 + a_block[5] 3484 1.00 + a_block[6] 2141 1.00 + a 1313 1.00 + bp 2439 1.00 + bpc 2609 1.00 + sigma_actor[1] 3874 1.00 + sigma_actor[2] 814 1.00 + sigma_actor[3] 951 1.00 + sigma_block[1] 1096 1.00 + sigma_block[2] 1356 1.00 + sigma_block[3] 1109 1.00 + Rho_actor[1,1] NaN NaN + Rho_actor[1,2] 3766 1.00 + Rho_actor[1,3] 3052 1.00 + Rho_actor[2,1] 3766 1.00 + Rho_actor[2,2] 10943 1.00 + Rho_actor[2,3] 4210 1.00 + Rho_actor[3,1] 3052 1.00 + Rho_actor[3,2] 4210 1.00 + Rho_actor[3,3] 8994 1.00 + Rho_block[1,1] NaN NaN + Rho_block[1,2] 3441 1.00 + Rho_block[1,3] 5484 1.00 + Rho_block[2,1] 3441 1.00 + Rho_block[2,2] 11102 1.00 + Rho_block[2,3] 4267 1.00 + Rho_block[3,1] 5484 1.00 + Rho_block[3,2] 4267 1.00 + Rho_block[3,3] 7014 1.00 + lp__ 372 1.01 +""" + diff --git a/scripts/non-identifiable.jl b/scripts/non-identifiable.jl new file mode 100644 index 00000000..2961f949 --- /dev/null +++ b/scripts/non-identifiable.jl @@ -0,0 +1,36 @@ +# ## Data + +using Distributions +using Random + +Random.seed!(1234) +y = rand(Normal(0,1), 100); + +# ## Model + +using Turing + +@model m8_4(y) = begin + ## Can't really set a Uniform[-Inf,Inf] on σ + α₁ ~ Uniform(-3000, 1000) + α₂ ~ Uniform(-1000, 3000) + σ ~ truncated(Cauchy(0,1), 0, Inf) + + y .~ Normal(α₁ + α₂, σ) +end + +chains = sample(m8_4(y), NUTS(0.65), 2000) + +# ## Output + +# \defaultoutput{} + +# ## Original output + +""" + mean sd 5.5% 94.5% n_eff Rhat + a1 -861.15 558.17 -1841.89 -31.04 7 1.43 + a2 861.26 558.17 31.31 1842.00 7 1.43 + sigma 0.97 0.07 0.89 1.09 9 1.17 +"""; + diff --git a/scripts/oceanic-tool-complexity.jl b/scripts/oceanic-tool-complexity.jl new file mode 100644 index 00000000..5ddebdd4 --- /dev/null +++ b/scripts/oceanic-tool-complexity.jl @@ -0,0 +1,56 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames + +delim = ';' +data_path = joinpath(TuringModels.project_root, "data", "Kline.csv") +df = CSV.read(data_path, DataFrame; delim) + +df.log_pop = log.(df.population) + +## New col contact_high, set binary values 1/0 if high/low contact +df.contact_high = [contact == "high" ? 1 : 0 for contact in df.contact] +df + +# ## Model + +using Turing + +## This is supposed to be a "bad" model since we take non-centered data for the +## predictor log_pop +@model m10_10stan(total_tools, log_pop, contact_high) = begin + α ~ Normal(0, 100) + βp ~ Normal(0, 1) + βc ~ Normal(0, 1) + βpc ~ Normal(0, 1) + + for i ∈ 1:length(total_tools) + λ = exp(α + βp*log_pop[i] + βc*contact_high[i] + + βpc*contact_high[i]*log_pop[i]) + total_tools[i] ~ Poisson(λ) + end +end; + +# ## Output + +chains = sample( + m10_10stan(df.total_tools, df.log_pop, df.contact_high), + NUTS(0.65), + 1000 +) + +# \defaultoutput{} + +# ## Original output + +m_10_10t_result = " + mean sd 5.5% 94.5% n_eff Rhat + a 0.94 0.37 0.36 1.53 3379 1 + bp 0.26 0.04 0.21 0.32 3347 1 + bc -0.08 0.84 -1.41 1.23 2950 1 + bpc 0.04 0.09 -0.10 0.19 2907 1 +"; + diff --git a/scripts/spatial-autocorrelation-oceanic.jl b/scripts/spatial-autocorrelation-oceanic.jl new file mode 100644 index 00000000..21362862 --- /dev/null +++ b/scripts/spatial-autocorrelation-oceanic.jl @@ -0,0 +1,73 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames +using RData +using LinearAlgebra + +data_dir = joinpath(TuringModels.project_root, "data") +kline_path = joinpath(data_dir, "Kline2.csv") +df = CSV.read(kline_path, DataFrame; delim=';') + +dmat_path = joinpath(data_dir, "islandsDistMatrix.rda") +dmat = load(dmat_path)["islandsDistMatrix"] + +df.society = 1:10 + +# ## Model + +using Turing + +@model function m13_7(Dmat, society, logpop, total_tools) + rhosq ~ truncated(Cauchy(0, 1), 0, Inf) + etasq ~ truncated(Cauchy(0, 1), 0, Inf) + bp ~ Normal(0, 1) + a ~ Normal(0, 10) + + ## GPL2 + SIGMA_Dmat = etasq * exp.(-rhosq * Dmat.^2) + SIGMA_Dmat = SIGMA_Dmat + 0.01I + SIGMA_Dmat = (SIGMA_Dmat' + SIGMA_Dmat) / 2 + g ~ MvNormal(zeros(10), SIGMA_Dmat) + + log_lambda = a .+ g[society] .+ bp * logpop + + total_tools .~ Poisson.(exp.(log_lambda)) +end + +chains = sample( + m13_7(dmat, df.society, df.logpop, df.total_tools), + NUTS(0.65), + 5000 +) + +# ## Output + +# \defaultoutput{} + +# ## Original output + +m_13_7_rethinking = """ +Inference for Stan model: 6422d8042e9cdd08dae2420ad26842f1. + 4 chains, each with iter=10000; warmup=2000; thin=1; + post-warmup draws per chain=8000, total post-warmup draws=32000. + + mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat + g[1] -0.27 0.01 0.44 -1.26 -0.50 -0.24 -0.01 0.53 3278 1 + g[2] -0.13 0.01 0.43 -1.07 -0.34 -0.10 0.12 0.67 3144 1 + g[3] -0.17 0.01 0.42 -1.10 -0.37 -0.14 0.07 0.59 3096 1 + g[4] 0.30 0.01 0.37 -0.51 0.12 0.30 0.50 1.02 3145 1 + g[5] 0.03 0.01 0.37 -0.77 -0.15 0.04 0.22 0.72 3055 1 + g[6] -0.46 0.01 0.38 -1.31 -0.64 -0.42 -0.23 0.19 3306 1 + g[7] 0.10 0.01 0.36 -0.70 -0.07 0.11 0.29 0.78 3175 1 + g[8] -0.26 0.01 0.37 -1.07 -0.43 -0.24 -0.06 0.40 3246 1 + g[9] 0.23 0.01 0.35 -0.52 0.07 0.25 0.42 0.89 3302 1 + g[10] -0.12 0.01 0.45 -1.04 -0.36 -0.11 0.13 0.77 5359 1 + a 1.31 0.02 1.15 -0.98 0.61 1.31 2.00 3.69 4389 1 + bp 0.25 0.00 0.11 0.02 0.18 0.24 0.31 0.47 5634 1 + etasq 0.34 0.01 0.53 0.03 0.10 0.20 0.39 1.52 5643 1 + rhosq 1.52 0.09 11.82 0.03 0.16 0.39 0.96 7.96 15955 1 + lp__ 925.98 0.03 2.96 919.16 924.20 926.34 928.14 930.67 7296 1 +"""; diff --git a/scripts/varying-intercepts-admission.jl b/scripts/varying-intercepts-admission.jl new file mode 100644 index 00000000..5a0be81a --- /dev/null +++ b/scripts/varying-intercepts-admission.jl @@ -0,0 +1,98 @@ +# ## Data + +import CSV +import TuringModels + +using DataFrames + +data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv") +df = CSV.read(data_path, DataFrame; delim=';') + +dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept))) +df.male = [g == "male" ? 1 : 0 for g in df.gender] +df.dept_id = [dept_map[de] for de in df.dept] + +# ## Model + +using Turing + +@model function m13_3(applications, dept_id, male, admit) + Rho ~ LKJ(2, 2.) + sigma_dept ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2) + bm ~ Normal(0, 1) + a ~ Normal(0, 1) + + dist_mu = [a, bm] + dist_Sigma = sigma_dept .* Rho .* sigma_dept' + dist_Sigma = (dist_Sigma' + dist_Sigma) / 2 + + a_bm_dept ~ filldist(MvNormal(dist_mu, dist_Sigma), 6) + + a_dept, bm_dept = a_bm_dept[1, :], a_bm_dept[2, :] + logit_p = a_dept[dept_id] + bm_dept[dept_id] .* male + + admit .~ BinomialLogit.(applications, logit_p) +end; + +# ## Output + +chains = sample( + m13_3(df.applications, df.dept_id, df.male, df.admit), + NUTS(0.95), + 1000 +) + +# \defaultoutput{} + +# ## Original output + +m_13_3_rethinking = """ +Inference for Stan model: f0d86ec689cbf7921aab4fc0f55616d2. + 4 chains, each with iter=5000; warmup=1000; thin=1; + post-warmup draws per chain=4000, total post-warmup draws=16000. + + mean se_mean sd 2.5% 25% 50% 75% + bm_dept[1] -0.79 0.00 0.27 -1.32 -0.97 -0.79 -0.61 + bm_dept[2] -0.21 0.00 0.33 -0.87 -0.42 -0.20 0.00 + bm_dept[3] 0.08 0.00 0.14 -0.18 -0.01 0.08 0.18 + bm_dept[4] -0.09 0.00 0.14 -0.37 -0.19 -0.09 0.00 + bm_dept[5] 0.12 0.00 0.19 -0.24 -0.01 0.12 0.25 + bm_dept[6] -0.12 0.00 0.27 -0.67 -0.29 -0.12 0.05 + a_dept[1] 1.30 0.00 0.25 0.82 1.13 1.30 1.47 + a_dept[2] 0.74 0.00 0.33 0.08 0.52 0.73 0.95 + a_dept[3] -0.65 0.00 0.09 -0.82 -0.71 -0.65 -0.59 + a_dept[4] -0.62 0.00 0.10 -0.83 -0.69 -0.62 -0.55 + a_dept[5] -1.13 0.00 0.11 -1.36 -1.21 -1.13 -1.05 + a_dept[6] -2.60 0.00 0.20 -3.01 -2.73 -2.60 -2.46 + a -0.49 0.01 0.73 -1.96 -0.91 -0.49 -0.07 + bm -0.16 0.00 0.24 -0.65 -0.29 -0.16 -0.02 + sigma_dept[1] 1.67 0.01 0.63 0.88 1.25 1.53 1.94 + sigma_dept[2] 0.50 0.00 0.26 0.16 0.33 0.45 0.61 + Rho[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 + Rho[1,2] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08 + Rho[2,1] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08 + Rho[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 + lp__ -2593.92 0.04 3.17 -2601.03 -2595.85 -2593.57 -2591.65 + 97.5% n_eff Rhat + bm_dept[1] -0.27 6918 1 + bm_dept[2] 0.46 9133 1 + bm_dept[3] 0.36 13016 1 + bm_dept[4] 0.19 12258 1 + bm_dept[5] 0.49 12953 1 + bm_dept[6] 0.40 12011 1 + a_dept[1] 1.81 6993 1 + a_dept[2] 1.39 9335 1 + a_dept[3] -0.48 13489 1 + a_dept[4] -0.42 12595 1 + a_dept[5] -0.91 14659 1 + a_dept[6] -2.23 13168 1 + a 0.96 9541 1 + bm 0.31 8946 1 + sigma_dept[1] 3.27 9058 1 + sigma_dept[2] 1.15 6626 1 + Rho[1,1] 1.00 NaN NaN + Rho[1,2] 0.44 11358 1 + Rho[2,1] 0.44 11358 1 + Rho[2,2] 1.00 16069 1 + lp__ -2588.70 5310 1 +"""; diff --git a/scripts/varying-intercepts-reedfrogs.jl b/scripts/varying-intercepts-reedfrogs.jl new file mode 100644 index 00000000..3ede73bd --- /dev/null +++ b/scripts/varying-intercepts-reedfrogs.jl @@ -0,0 +1,84 @@ +# ## Data + +import CSV + +using DataFrames +using TuringModels + +data_path = joinpath(TuringModels.project_root, "data", "reedfrogs.csv") +df = CSV.read(data_path, DataFrame; delim=';'); +@assert size(df) == (48, 5) ## hide +df.tank = 1:nrow(df) + +# ## Model + +using Turing + +@model function reedfrogs(density, tank, surv, n) + a_tank ~ filldist(Normal(0, 1.5), n) + + logitp = a_tank[tank] + surv .~ BinomialLogit.(density, logitp) +end; + +# ## Output + +n = nrow(df) +model = reedfrogs(df.density, df.tank, df.surv, n) +chains = sample(model, NUTS(0.65), 1000) + +# \defaultoutput{} + +# ## Original output (Edition 1) + +""" + mean sd 5.5% 94.5% n_eff Rhat +a_tank[1] 2.49 1.16 0.85 4.53 1079 1 +a_tank[2] 5.69 2.75 2.22 10.89 1055 1 +a_tank[3] 0.89 0.75 -0.23 2.16 1891 1 +a_tank[4] 5.71 2.70 2.21 10.85 684 1 +a_tank[5] 2.52 1.14 0.92 4.42 1640 1 +a_tank[6] 2.49 1.13 0.94 4.52 1164 1 +a_tank[7] 5.74 2.71 2.25 10.86 777 1 +a_tank[8] 2.52 1.19 0.95 4.42 1000 1 +a_tank[9] -0.46 0.69 -1.62 0.55 2673 1 +a_tank[10] 2.53 1.19 0.93 4.59 1430 1 +a_tank[11] 0.93 0.72 -0.17 2.11 1387 1 +a_tank[12] 0.47 0.74 -0.63 1.70 1346 1 +a_tank[13] 0.91 0.76 -0.25 2.30 1559 1 +a_tank[14] 0.00 0.66 -1.04 1.06 2085 1 +a_tank[15] 2.50 1.19 0.95 4.40 1317 1 +a_tank[16] 2.50 1.14 0.98 4.31 1412 1 +a_tank[17] 3.49 1.12 1.94 5.49 945 1 +a_tank[18] 2.59 0.75 1.50 3.81 1561 1 +a_tank[19] 2.11 0.64 1.15 3.15 1712 1 +a_tank[20] 6.40 2.57 3.11 11.04 996 1 +a_tank[21] 2.59 0.74 1.54 3.93 1233 1 +a_tank[22] 2.63 0.79 1.49 4.01 1184 1 +a_tank[23] 2.64 0.83 1.45 4.13 1379 1 +a_tank[24] 1.74 0.59 0.85 2.72 1736 1 +a_tank[25] -1.19 0.45 -1.90 -0.50 2145 1 +a_tank[26] 0.09 0.41 -0.53 0.78 2167 1 +a_tank[27] -1.75 0.56 -2.65 -0.88 1666 1 +a_tank[28] -0.58 0.43 -1.25 0.08 1567 1 +a_tank[29] 0.08 0.39 -0.54 0.71 3053 1 +a_tank[30] 1.43 0.49 0.66 2.24 2754 1 +a_tank[31] -0.79 0.44 -1.50 -0.12 1299 1 +a_tank[32] -0.42 0.41 -1.12 0.23 1661 1 +a_tank[33] 3.84 1.08 2.31 5.70 808 1 +a_tank[34] 3.00 0.85 1.83 4.36 1038 1 +a_tank[35] 2.96 0.82 1.82 4.25 1578 1 +a_tank[36] 2.14 0.55 1.31 3.08 1734 1 +a_tank[37] 2.12 0.56 1.31 3.04 1131 1 +a_tank[38] 6.72 2.62 3.45 11.44 706 1 +a_tank[39] 2.95 0.73 1.85 4.08 1509 1 +a_tank[40] 2.48 0.65 1.53 3.61 1731 1 +a_tank[41] -2.15 0.57 -3.11 -1.29 1231 1 +a_tank[42] -0.67 0.35 -1.22 -0.14 1444 1 +a_tank[43] -0.54 0.35 -1.12 0.03 1776 1 +a_tank[44] -0.43 0.34 -1.00 0.10 1735 1 +a_tank[45] 0.54 0.36 -0.04 1.14 1376 1 +a_tank[46] -0.67 0.34 -1.25 -0.15 1619 1 +a_tank[47] 2.14 0.55 1.31 3.04 1916 1 +a_tank[48] -0.06 0.35 -0.61 0.50 1932 1 +"""; diff --git a/scripts/varying-slopes-cafe.jl b/scripts/varying-slopes-cafe.jl new file mode 100644 index 00000000..087c8450 --- /dev/null +++ b/scripts/varying-slopes-cafe.jl @@ -0,0 +1,109 @@ +# ## Data + +import CSV + +using DataFrames +using TuringModels + +data_path = joinpath(TuringModels.project_root, "data", "d_13_1.csv") +df = CSV.read(data_path, DataFrame); + +## DataFrame `df` is shown Section [df](#df). + +# ## Model + +using Turing + +@model function m13_1(cafe, afternoon, wait) + Rho ~ LKJ(2, 1.) + sigma ~ truncated(Cauchy(0, 2), 0, Inf) + sigma_cafe ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2) + a ~ Normal(0, 10) + b ~ Normal(0, 10) + + dist_mu = [a, b] + dist_Sigma = sigma_cafe .* Rho .* sigma_cafe' + dist_Sigma = (dist_Sigma' + dist_Sigma) / 2 + a_b_cafe ~ filldist(MvNormal(dist_mu, dist_Sigma), 20) + + a_cafe = a_b_cafe[1, :] + b_cafe = a_b_cafe[2, :] + + μ = a_cafe[cafe] + b_cafe[cafe] .* afternoon + wait .~ Normal.(μ, sigma) +end; + +# ## Output + +chains = sample( + m13_1(df.cafe, df.afternoon, df.wait), + ## This model fails on NUTS(0.65). + Turing.NUTS(0.95), + 1000 +) + +# \defaultoutput{} + +# ## Original output + +""" +Inference for Stan model: a73b0bd01032773825c6abf5575fd6e4. + 2 chains, each with iter=5000; warmup=2000; thin=1; + post-warmup draws per chain=3000, total post-warmup draws=6000. + + mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat + b_cafe[1] -1.29 0.00 0.18 -1.69 -1.40 -1.28 -1.18 -0.96 2548 1.00 + b_cafe[2] -1.20 0.00 0.18 -1.57 -1.31 -1.21 -1.09 -0.83 3288 1.00 + b_cafe[3] -1.26 0.00 0.18 -1.63 -1.36 -1.25 -1.15 -0.91 4302 1.00 + b_cafe[4] -1.29 0.00 0.18 -1.68 -1.39 -1.28 -1.18 -0.96 2960 1.00 + b_cafe[5] -1.26 0.00 0.20 -1.68 -1.38 -1.25 -1.14 -0.88 3406 1.00 + b_cafe[6] -1.28 0.00 0.18 -1.66 -1.38 -1.27 -1.17 -0.93 2983 1.00 + b_cafe[7] -1.23 0.00 0.18 -1.62 -1.33 -1.23 -1.13 -0.86 4430 1.00 + b_cafe[8] -1.25 0.00 0.17 -1.62 -1.35 -1.24 -1.15 -0.90 3838 1.00 + b_cafe[9] -1.13 0.01 0.19 -1.45 -1.26 -1.16 -1.02 -0.70 1372 1.00 + b_cafe[10] -1.19 0.00 0.18 -1.54 -1.30 -1.20 -1.09 -0.80 3924 1.00 + b_cafe[11] -1.04 0.01 0.22 -1.38 -1.20 -1.07 -0.90 -0.55 416 1.01 + b_cafe[12] -1.21 0.00 0.18 -1.55 -1.32 -1.22 -1.11 -0.82 3781 1.00 + b_cafe[13] -1.32 0.00 0.19 -1.76 -1.43 -1.30 -1.20 -0.99 1509 1.00 + b_cafe[14] -1.37 0.01 0.20 -1.82 -1.49 -1.34 -1.23 -1.04 760 1.01 + b_cafe[15] -1.52 0.02 0.27 -2.11 -1.70 -1.49 -1.31 -1.13 161 1.02 + b_cafe[16] -1.18 0.00 0.18 -1.51 -1.29 -1.20 -1.08 -0.79 3430 1.00 + b_cafe[17] -1.16 0.00 0.19 -1.50 -1.28 -1.18 -1.05 -0.73 2034 1.00 + b_cafe[18] -1.28 0.00 0.21 -1.72 -1.41 -1.28 -1.16 -0.86 3166 1.00 + b_cafe[19] -1.05 0.01 0.22 -1.38 -1.21 -1.09 -0.91 -0.55 363 1.01 + b_cafe[20] -1.10 0.01 0.20 -1.43 -1.23 -1.12 -0.98 -0.64 613 1.01 + a_cafe[1] 4.08 0.00 0.18 3.74 3.97 4.08 4.20 4.44 5407 1.00 + a_cafe[2] 2.38 0.00 0.18 2.03 2.26 2.38 2.50 2.72 5091 1.00 + a_cafe[3] 3.94 0.00 0.18 3.60 3.82 3.94 4.06 4.30 6412 1.00 + a_cafe[4] 3.45 0.00 0.18 3.11 3.34 3.45 3.57 3.81 5699 1.00 + a_cafe[5] 2.14 0.00 0.18 1.79 2.02 2.14 2.26 2.50 5380 1.00 + a_cafe[6] 4.26 0.00 0.17 3.92 4.15 4.26 4.38 4.61 5192 1.00 + a_cafe[7] 3.56 0.00 0.18 3.21 3.44 3.56 3.68 3.91 5495 1.00 + a_cafe[8] 3.79 0.00 0.18 3.44 3.68 3.79 3.91 4.14 5661 1.00 + a_cafe[9] 3.89 0.00 0.18 3.53 3.77 3.89 4.01 4.23 4135 1.00 + a_cafe[10] 3.69 0.00 0.18 3.34 3.57 3.69 3.81 4.05 5761 1.00 + a_cafe[11] 2.47 0.00 0.19 2.11 2.35 2.48 2.60 2.84 1548 1.00 + a_cafe[12] 4.08 0.00 0.17 3.74 3.97 4.08 4.20 4.42 5712 1.00 + a_cafe[13] 3.88 0.00 0.18 3.52 3.75 3.87 4.00 4.24 4061 1.00 + a_cafe[14] 3.33 0.00 0.19 2.97 3.21 3.33 3.46 3.70 3082 1.00 + a_cafe[15] 4.23 0.01 0.20 3.85 4.09 4.22 4.36 4.65 471 1.01 + a_cafe[16] 3.60 0.00 0.17 3.25 3.48 3.60 3.71 3.93 5527 1.00 + a_cafe[17] 4.43 0.00 0.18 4.09 4.31 4.43 4.55 4.78 5089 1.00 + a_cafe[18] 6.10 0.00 0.19 5.73 5.97 6.09 6.22 6.46 4780 1.00 + a_cafe[19] 3.50 0.00 0.19 3.12 3.38 3.50 3.63 3.86 1800 1.00 + a_cafe[20] 3.90 0.00 0.18 3.53 3.78 3.90 4.03 4.25 3465 1.00 + a 3.73 0.00 0.21 3.32 3.60 3.73 3.87 4.16 7131 1.00 + b -1.23 0.00 0.09 -1.40 -1.29 -1.23 -1.18 -1.06 2021 1.00 + sigma_cafe[1] 0.91 0.00 0.17 0.65 0.80 0.89 1.01 1.29 5579 1.00 + sigma_cafe[2] 0.21 0.01 0.12 0.01 0.12 0.20 0.29 0.46 72 1.05 + sigma 0.49 0.00 0.03 0.44 0.47 0.49 0.51 0.55 2271 1.00 + Rho[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN + Rho[1,2] -0.17 0.01 0.34 -0.75 -0.42 -0.20 0.06 0.57 3300 1.00 + Rho[2,1] -0.17 0.01 0.34 -0.75 -0.42 -0.20 0.06 0.57 3300 1.00 + Rho[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 5777 1.00 + lp__ 62.41 3.30 17.59 41.97 52.03 58.04 66.26 118.46 28 1.12 +"""; + +## df + +df diff --git a/scripts/weakly-informative-priors.jl b/scripts/weakly-informative-priors.jl new file mode 100644 index 00000000..f5f8a7fa --- /dev/null +++ b/scripts/weakly-informative-priors.jl @@ -0,0 +1,29 @@ +# ## Data + +y = [-1,1] + +# ## Model + +using Turing + +@model function m8_3(y) + α ~ Normal(1, 10) + σ ~ truncated(Cauchy(0, 1), 0, Inf) + + μ = α + y .~ Normal(μ, σ) +end; + +# ## Output + +chains = sample(m8_3(y), NUTS(0.65), 1000) + +# \defaultoutput{} + +# ## Original output + +""" + mean sd 5.5% 94.5% n_eff Rhat +alpha 0.09 1.63 -2.13 2.39 959 1 +sigma 2.04 2.05 0.68 4.83 1090 1 +"""; diff --git a/scripts/wild-chain.jl b/scripts/wild-chain.jl new file mode 100644 index 00000000..ffbbc836 --- /dev/null +++ b/scripts/wild-chain.jl @@ -0,0 +1,21 @@ +# ## Data + +y = [-1, 1] + +# ## Model + +using Turing + +@model function m8_2(y) + σ ~ FlatPos(0.0) ## improper prior with probability one everywhere above 0.0 + α ~ Flat() ## improper prior with pobability one everywhere + + y .~ Normal(α, σ) +end + +# ## Output + +chains = sample(m8_2(y), NUTS(0.65), 1000) + +# \defaultoutput{} +