Skip to content
This repository has been archived by the owner on Sep 10, 2023. It is now read-only.

Commit

Permalink
Move more models to Literate (#58)
Browse files Browse the repository at this point in the history
* Add compat entry for StatsFuns

* Let height correspond to Edition 2

* Fix ignoring-gender-admit

* Fix Africa
  • Loading branch information
rikhuijzer committed Jan 5, 2021
1 parent 18c3604 commit a20059d
Show file tree
Hide file tree
Showing 15 changed files with 274 additions and 282 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ CSV = "0.8"
DataFrames = "0.22"
Franklin = "0.10"
NodeJS = "1.1"
StatsFuns = "0.9"
StatsPlots = "0.14"
Turing = "0.15"

8 changes: 8 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,14 @@ This package contains Julia versions of the mcmc models contained in the R packa

This package implements the models using [TuringLang/Turing.jl](https://github.com/TuringLang).

## Usage

Most of the scripts and output can be inspected via the webpages.
If you want to run the scripts yourselves, then you can either

1. copy the code from the webpages and the data from this repository, and run the scripts **or**
1. clone this repository and run one of the files in `scripts`. For example, `julia -i scripts/basic-example.jl`.

## Versions

### v1.1.2
Expand Down
2 changes: 1 addition & 1 deletion _layout/pgwrap.html
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,6 @@
<div class="main-header">
<a id="github"
target="_blank"
href="//github.com/StatisticalRethinkingJulia/TuringModels.jl">TuringModels on GitHub
href="https://github.com/StatisticalRethinkingJulia/TuringModels.jl">TuringModels on GitHub
</a>
</div>
4 changes: 2 additions & 2 deletions config.md
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
@def prepath = "TuringModels.jl"
@def website_title = "TuringModels"
@def website_descr = "Statistical Rethinking models in Julia"
@def website_url = "https://rikhuijzer.github.io/TuringModels.jl/"
@def website_url = "https://statisticalrethinkingjulia.github.io/TuringModels.jl/"

@def author = "Rob J Goedman and contributors"
@def author = "Rob J Goedman, Rik Huijzer and contributors"

@def mintoclevel = 2

Expand Down
4 changes: 3 additions & 1 deletion index.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,10 @@ The code for this website is available on [GitHub](https://github.com/Statistica

Each page aims to contain all the code required to reproduce the results.
In other words, it should be possible to get the same output by just copying the code and the accompanying dataset.
Furthermore, we use unicode symbols in the models.
Furthermore, we try to stick to Julia styling conventions where possible.
Therefore, we use unicode symbols in the models.
For example, where the book lists `alpha`, we will use `α` because the Julia language allows that.
As another example, we call DataFrame variables `df`.

### Version

Expand Down
30 changes: 15 additions & 15 deletions models.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,41 +5,41 @@ rss = "A list of the Turing versions of the Statistical Rethinking models."
tags = ["models", "statistics"]
+++

This page lists the available Turing models for different versions of the Statistical Rethinking book.
Some pages list code, tables and plots, and others show only code.

This page lists the available Turing models for different versions of the Statistical Rethinking books.

\toc

Before you look at the models below, you might want to look at the [Basic Example](basic-example).

## 2nd Edition (2020)

- [globe.qa](globe-tossing): Globe tossing
- [m4.1](height): Gaussian model of height
- [m8.1](africa): Africa
- [m12.1](beta-binomial): Beta-binomial
- [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]

## 1st Edition (2015)

- [m2.1](globe-tossing): Globe tossing
- [m4.1](height): Gaussian model of height
- [m8.1](africa): Africa
- [m8.2](wild-chain): Wild chain
- [m8.3](weakly-informative-priors): Weakly informative priors
- [m8.4](non-identifiable): Non-identifiable model [broken]
- [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]
- [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 [not-implemented]
- [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]
- [m14.1](varying-slopes-cafe): Varying slopes cafes

## 2nd Edition (2020)

- [globe.qa](globe-tossing): Globe tossing
- [m4.2](height): Gaussian model of height
- [m8.1](africa): Africa
- [m12.1](beta-binomial): Beta-binomial
- [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]
91 changes: 3 additions & 88 deletions models/africa.md
Original file line number Diff line number Diff line change
@@ -1,89 +1,4 @@
+++
title = "Africa"
+++
@def title = "Africa"
@def showall = true

This is the first Stan model in Statistical Rethinking Edition 1 (page 249).

\toc

## Data

```julia:data
using DataFrames
using Distributions
using Turing
using TuringModels # hide
import CSV
data_path = joinpath(TuringModels.project_root, "data", "rugged.csv")
df = CSV.read(data_path, DataFrame)
df.log_gdp = [log(x) for x in df.rgdppc_2000]
TODO = "Replace by dropmissing or something similar."
notisnan(e) = !ismissing(e)
df = df[map(notisnan, df[:, :rgdppc_2000]), :];
df = select(df, :log_gdp, :rugged, :cont_africa)
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
@model function model_fn(y, x₁, x₂)
σ ~ truncated(Cauchy(0, 2), 0, Inf)
βAR ~ Normal(0, 10)
βR ~ Normal(0, 10)
βA ~ Normal(0, 10)
α ~ Normal(0, 100)
μ = α .+ βR * x₁ .+ βA * x₂ .+ βAR * x₁ .* x₂
y .~ Normal.(μ, σ)
end
model = model_fn(df.log_gdp, df.rugged, df.cont_africa)
chains = sample(model, NUTS(0.65), 1000)
```
\output{model}

## Output

\defaultoutput{}

## Original output

```
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 9.22360053 0.139119116 0.0021996664 0.0034632816 1000
bR -0.20196346 0.076106388 0.0012033477 0.0018370185 1000
bA -1.94430980 0.227080488 0.0035904578 0.0057840746 1000
bAR 0.39071684 0.131889143 0.0020853505 0.0032749642 1000
sigma 0.95036370 0.052161768 0.0008247500 0.0009204073 1000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
a 8.95307475 9.12719750 9.2237750 9.31974000 9.490234250
bR -0.35217930 -0.25334425 -0.2012855 -0.15124725 -0.054216855
bA -2.39010825 -2.09894500 -1.9432550 -1.78643000 -1.513974250
bAR 0.13496995 0.30095575 0.3916590 0.47887625 0.650244475
sigma 0.85376115 0.91363250 0.9484920 0.98405750 1.058573750
```

## Appendix

### Appendix I

\tableinput{}{./code/output/data.csv}
\literate{/scripts/africa.jl}
32 changes: 2 additions & 30 deletions models/basic-example.md
Original file line number Diff line number Diff line change
@@ -1,34 +1,6 @@
+++
title = "Basic Example"
showall = true
+++

\toc

## Model definition

We define a simple Normal model with unknown mean and variance.

```julia:model
import CSV
using DataFrames
using StatsPlots
using Turing
@model gdemo(x, y) = begin
s ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s))
x ~ Normal(m, sqrt(s))
y ~ Normal(m, sqrt(s))
end
```

and run the sampler:

```julia:sampler
chains = sample(gdemo(1.5, 2), NUTS(0.65), 1000)
```

## Output

\defaultoutput{}
\literate{/scripts/basic-example.jl}
1 change: 0 additions & 1 deletion models/globe-tossing.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
+++
title = "Globe tossing"
showall = true
reeval = true
+++

\literate{/scripts/globe-tossing.jl}
96 changes: 11 additions & 85 deletions models/height.md
Original file line number Diff line number Diff line change
@@ -1,103 +1,29 @@
+++
title = "Gaussian model of height"
showall = true
+++

This model is based on data of the !Kung San people \citep{pHowell2010}.
The prior assumes that the height of a man is 178 cm because McElreath has that lenght, so it is a reasonable prior for the height.
The prior assumes that the height of a man is most likely to be 178 cm because McElreath has that lenght.

McElreath defines the model as
In Edition 2, McElreath defines the model as

$$
\begin{aligned}
h_i &\sim \text{Normal}(\mu, \sigma) \\
\mu &\sim \text{Normal}(178, 20) \\
\sigma &\sim \text{Uniform}(0, 50)
\end{aligned}
\begin{aligned}
h_i &\sim \text{Normal}(\mu, \sigma) \\
\mu &\sim \text{Normal}(178, 0.1) \\
\sigma &\sim \text{Uniform}(0, 50)
\end{aligned}
$$

\toc
\toc

## Model definition

```julia:height
import CSV
using DataFrames
using StatsPlots
using Turing
using TuringModels
data_path = joinpath(TuringModels.project_root, "data", "Howell1.csv")
df = CSV.read(data_path, DataFrame; delim=';')
# Use only adults and center the weight observations
df2 = filter(row -> row.age >= 18, df)
mean_weight = mean(df2.weight)
df2.weight_c = df2.weight .- mean_weight
@model function line(x, y)
alpha ~ Normal(178.0, 100.0)
beta ~ Normal(0.0, 10.0)
sigma ~ Uniform(0, 50)
mu = alpha .+ beta*x
y .~ Normal.(mu, sigma)
end
x = df2.weight_c
y = df2.height
chains = sample(line(x, y), NUTS(0.65), 1000)
```
\output{height}

```julia:write_helper
# hideall
output_dir = @OUTPUT
function write_svg(name, p)
fig_path = joinpath(output_dir, "$name.svg")
StatsPlots.savefig(fig_path)
end
```
\output{write_helper}

## Output

\defaultoutput{}

## Original output

```
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
alpha 154.597086 0.27326431 0.0043206882 0.0036304132 1000
beta 0.906380 0.04143488 0.0006551430 0.0006994720 1000
sigma 5.106643 0.19345409 0.0030587777 0.0032035103 1000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
alpha 154.0610000 154.4150000 154.5980000 154.7812500 155.1260000
beta 0.8255494 0.8790695 0.9057435 0.9336445 0.9882981
sigma 4.7524368 4.9683400 5.0994450 5.2353100 5.5090128
```
\literate{/scripts/height.jl}

## References

\biblabel{pHowell2010}{Howell, 2010}
Howell, N. (2010).
Life Histories of the Dobe !Kung: Food, Fatness, and Well-being over the Life-span.
Origins of Human Behavior and Culture.
University of California Press

## Appendix
### df2
```julia:data
write_csv(name, data) = CSV.write(joinpath(@OUTPUT, "$name.csv"), data) # hide
write_csv("data", # hide
df2
) # hide
```
\output{data}
\tableinput{}{./code/output/data.csv}
Loading

0 comments on commit a20059d

Please sign in to comment.