From 30e59f04e416bb18abe324831f62dec3b983d066 Mon Sep 17 00:00:00 2001 From: Penelope Yong Date: Fri, 30 Aug 2024 23:15:09 +0100 Subject: [PATCH 1/3] Expand 'Getting Started' page - Bump the minimum version of Julia. - Explain the significance of parameters to @model function. - Explain how to obtain the actual samples in the chain. - Expand on the maths involved in the analytical posterior calculation (otherwise, to someone who hasn't seen it before, the formulas look too magic). --- tutorials/docs-00-getting-started/index.qmd | 119 ++++++++++++++------ 1 file changed, 86 insertions(+), 33 deletions(-) diff --git a/tutorials/docs-00-getting-started/index.qmd b/tutorials/docs-00-getting-started/index.qmd index 17197f977..f2d949345 100644 --- a/tutorials/docs-00-getting-started/index.qmd +++ b/tutorials/docs-00-getting-started/index.qmd @@ -16,13 +16,9 @@ Pkg.instantiate(); To use Turing, you need to install Julia first and then install Turing. -### Install Julia +You will need to install Julia 1.7 or greater, which you can get from [the official Julia website](http://julialang.org/downloads/). -You will need to install Julia 1.3 or greater, which you can get from [the official Julia website](http://julialang.org/downloads/). - -### Install Turing.jl - -Turing is an officially registered Julia package, so you can install a stable version of Turing by running the following in the Julia REPL: +Turing is officially registered in the [Julia General package registry](https://github.com/JuliaRegistries/General), which means that you can install a stable version of Turing by running the following in the Julia REPL: ```{julia} #| output: false @@ -30,82 +26,139 @@ using Pkg Pkg.add("Turing") ``` -You can check if all tests pass by running `Pkg.test("Turing")` (it might take a long time) +### Example usage -### Example +Here is a simple example showing Turing in action. -Here's a simple example showing Turing in action. - -First, we can load the Turing and StatsPlots modules +First, we load the Turing and StatsPlots modules. +The latter is required for visualising the results. ```{julia} using Turing using StatsPlots ``` -Then, we define a simple Normal model with unknown mean and variance +Then, we define a simple Normal model with unknown mean and variance. +Here, both `x` and `y` are observed values (since they are function parameters), whereas `m` and `s²` are the parameters to be inferred. ```{julia} @model function gdemo(x, y) s² ~ InverseGamma(2, 3) m ~ Normal(0, sqrt(s²)) x ~ Normal(m, sqrt(s²)) - return y ~ Normal(m, sqrt(s²)) + y ~ Normal(m, sqrt(s²)) end ``` -Then we can run a sampler to collect results. In this case, it is a Hamiltonian Monte Carlo sampler +Suppose we have some data `x` and `y` that we want to infer the mean and variance for. +We can pass these data as arguments to the `gdemo` function, and run a sampler to collect the results. +In this specific example, we collect 1000 samples using the No U-Turn Sampler (NUTS) algorithm, which is a Hamiltonian Monte Carlo method. ```{julia} chn = sample(gdemo(1.5, 2), NUTS(), 1000, progress=false) ``` -We can plot the results +We can plot the results: ```{julia} plot(chn) ``` -In this case, because we use the normal-inverse gamma distribution as a conjugate prior, we can compute its updated mean as follows: +and obtain summary statistics by indexing the chain: + +```{julia} +mean(chn[:m]), mean(chn[:s²]) +``` + + +### Verifying the results + +To verify the results of this example, we can analytically solve for the posterior distribution. +Our prior is a normal-inverse gamma distribution: + +$$\begin{align} + P(s^2) &= \Gamma^{-1}(s^2 | \alpha = 2, \beta = 3) \\ + P(m) &= \mathcal{N}(m | \mu = 0, \sigma^2 = s^2) \\ +\implies P(m, s^2) &= \text{N-}\Gamma^{-1}(m, s^2 | \mu = 0, \lambda = 1, \alpha = 2, \beta = 3). +\end{align}$$ + +This is a conjugate prior for a normal distribution with unknown mean and variance. +So, the posterior distribution given some data $\mathbf{x}$ is also a normal-inverse gamma distribution, which we can compute analytically: + +$$P(m, s^2 | \mathbf{x}) = \text{N-}\Gamma^{-1}(m, s^2 | \mu', \lambda', \alpha', \beta').$$ + +The four updated parameters are given respectively (see for example [Wikipedia](https://en.wikipedia.org/wiki/Normal-gamma_distribution#Posterior_distribution_of_the_parameters)) by + +$$\begin{align} +\mu' &= \frac{\lambda\mu + n\bar{x}}{\lambda + n}; & +\lambda' &= \lambda + n; \\[0.3cm] +\alpha' &= \alpha + \frac{n}{2}; & +\beta' &= \beta + \frac{1}{2}\left(nv + \frac{\lambda n(\bar{x} - \mu)^2}{\lambda + n}\right), +\end{align}$$ + +where $n$ is the number of data points, $\bar{x}$ the mean of the data, and $v$ the uncorrected sample variance of the data. ```{julia} s² = InverseGamma(2, 3) -m = Normal(0, 1) +α, β = shape(s²), scale(s²) + +# We defined m ~ Normal(0, sqrt(s²)), and the normal-inverse gamma +# distribution is defined by m ~ Normal(µ, sqrt(s²/λ)), which lets +# us identify the following +µ, λ = 0, 1 + data = [1.5, 2] x_bar = mean(data) -N = length(data) +n = length(data) +v = var(data; corrected=false) + +µ´ = ((λ * µ) + (n * x_bar)) / (λ + n) +λ´ = λ + n +α´ = α + n / 2 +β´ = β + ((n * v) + ((λ * n * (x_bar - µ)^2) / (λ + n))) / 2 -mean_exp = (m.σ * m.μ + N * x_bar) / (m.σ + N) +(µ´, λ´, α´, β´) ``` -We can also compute the updated variance +We can use these to analytically calculate, for example, the expectation values of the mean and variance parameters in the posterior distribution: $E[m] = \mu'$, and $E[s^2] = \beta'/(\alpha + - 1)$. ```{julia} -updated_alpha = shape(s²) + (N / 2) -updated_beta = - scale(s²) + - (1 / 2) * sum((data[n] - x_bar)^2 for n in 1:N) + - (N * m.σ) / (N + m.σ) * ((x_bar)^2) / 2 -variance_exp = updated_beta / (updated_alpha - 1) +mean_exp = µ´ +variance_exp = β´ / (α´ - 1) + +(mean_exp, variance_exp) ``` -Finally, we can check if these expectations align with our HMC approximations from earlier. We can compute samples from a normal-inverse gamma following the equations given [here](https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution#Generating_normal-inverse-gamma_random_variates). +Alternatively, we can compare the two distributions visually (one parameter at a time). +To do so, we will directly sample from the analytical posterior distribution using the procedure described in [Wikipedia](https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution#Generating_normal-inverse-gamma_random_variates). ```{julia} -function sample_posterior(alpha, beta, mean, lambda, iterations) +function sample_posterior_analytic(µ´, λ´, α´, β´, iterations) samples = [] for i in 1:iterations - sample_variance = rand(InverseGamma(alpha, beta), 1) - sample_x = rand(Normal(mean, sqrt(sample_variance[1]) / lambda), 1) - samples = append!(samples, sample_x) + sampled_s² = rand(InverseGamma(α´, β')) + sampled_m = rand(Normal(µ´, sqrt(sampled_s² / λ´))) + samples = push!(samples, (sampled_m, sampled_s²)) end return samples end -analytical_samples = sample_posterior(updated_alpha, updated_beta, mean_exp, 2, 1000); +analytical_samples = sample_posterior_analytic(µ´, λ´, α´, β´, 1000); ``` +Comparing the distribution of the mean: + ```{julia} -density(analytical_samples; label="Posterior (Analytical)") +density([s[1] for s in analytical_samples]; label="Posterior (Analytical)") density!(chn[:m]; label="Posterior (HMC)") ``` + +and the distribution of the variance: + +```{julia} +density([s[2] for s in analytical_samples]; label="Posterior (Analytical)") +density!(chn[:s²]; label="Posterior (HMC)") +``` + +we see that our HMC sampler has given us a good representation of the posterior distribution. From 257c8c6ae48c40384aae3421c65aa58bde249d62 Mon Sep 17 00:00:00 2001 From: Penelope Yong Date: Tue, 10 Sep 2024 16:33:39 +0100 Subject: [PATCH 2/3] Remove quick start -- it's repeated content --- _quarto.yml | 16 ++-- .../docs-12-using-turing-guide/index.qmd | 4 +- .../index.qmd | 74 ------------------- 3 files changed, 11 insertions(+), 83 deletions(-) delete mode 100755 tutorials/docs-14-using-turing-quick-start/index.qmd diff --git a/_quarto.yml b/_quarto.yml index bbc4c44f1..52d1d2927 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -50,15 +50,14 @@ website: - text: documentation collapse-level: 1 contents: - - section: "Documentation" + - section: "For Users" # href: tutorials/index.qmd, This page will be added later so keep this line commented contents: - - section: "Using Turing - Modelling Syntax and Interface" + - section: "Using the Turing library" collapse-level: 1 contents: - tutorials/docs-00-getting-started/index.qmd - - text: "Quick Start" - href: tutorials/docs-14-using-turing-quick-start/index.qmd + - tutorials/00-introduction/index.qmd - tutorials/docs-12-using-turing-guide/index.qmd - text: "Mode Estimation" href: tutorials/docs-17-mode-estimation/index.qmd @@ -70,9 +69,8 @@ website: - text: "External Samplers" href: tutorials/docs-16-using-turing-external-samplers/index.qmd - - section: "Using Turing - Tutorials" + - section: "Examples of Turing Models" contents: - - tutorials/00-introduction/index.qmd - text: Gaussian Mixture Models href: tutorials/01-gaussian-mixture-model/index.qmd - tutorials/02-logistic-regression/index.qmd @@ -97,13 +95,15 @@ website: - text: "Gaussian Process Latent Variable Models" href: tutorials/12-gplvm/index.qmd - - section: "Developers: Contributing" + - section: "For Developers" + contents: + - section: "Contributing" collapse-level: 1 contents: - text: "How to Contribute" href: tutorials/docs-01-contributing-guide/index.qmd - - section: "Developers: PPL" + - section: "How Turing Works" collapse-level: 1 contents: - tutorials/docs-05-for-developers-compiler/index.qmd diff --git a/tutorials/docs-12-using-turing-guide/index.qmd b/tutorials/docs-12-using-turing-guide/index.qmd index b66b4d29d..56d59a144 100755 --- a/tutorials/docs-12-using-turing-guide/index.qmd +++ b/tutorials/docs-12-using-turing-guide/index.qmd @@ -1,5 +1,5 @@ --- -title: Guide +title: "Turing's Core Functionality" engine: julia --- @@ -10,6 +10,8 @@ using Pkg; Pkg.instantiate(); ``` +This article provides an overview of the core functionality in Turing.jl, which are likely to be used across a wide range of models. + ## Basics ### Introduction diff --git a/tutorials/docs-14-using-turing-quick-start/index.qmd b/tutorials/docs-14-using-turing-quick-start/index.qmd deleted file mode 100755 index 67d848752..000000000 --- a/tutorials/docs-14-using-turing-quick-start/index.qmd +++ /dev/null @@ -1,74 +0,0 @@ ---- -pagetitle: Quick Start -title: Probabilistic Programming in Thirty Seconds -engine: julia ---- - -```{julia} -#| echo: false -#| output: false -using Pkg; -Pkg.instantiate(); -``` - -If you are already well-versed in probabilistic programming and want to take a quick look at how Turing's syntax works or otherwise just want a model to start with, we have provided a complete Bayesian coin-flipping model below. - -This example can be run wherever you have Julia installed (see [Getting Started](../docs-00-getting-started/index.qmd), but you will need to install the packages `Turing` and `StatsPlots` if you have not done so already. - -This is an excerpt from a more formal example which can be found [here](../00-introduction/index.qmd). - -## Import Libraries -```{julia} -# Import libraries. -using Turing, StatsPlots, Random -``` - -```{julia} -# Set the true probability of heads in a coin. -p_true = 0.5 - -# Iterate from having seen 0 observations to 100 observations. -Ns = 0:100 -``` - -```{julia} -# Draw data from a Bernoulli distribution, i.e. draw heads or tails. -Random.seed!(12) -data = rand(Bernoulli(p_true), last(Ns)) -``` - - -## Declare Turing Model -```{julia} -# Declare our Turing model. -@model function coinflip(y) - # Our prior belief about the probability of heads in a coin. - p ~ Beta(1, 1) - - # The number of observations. - N = length(y) - for n in 1:N - # Heads or tails of a coin are drawn from a Bernoulli distribution. - y[n] ~ Bernoulli(p) - end -end -``` - - -## Setting HMC Sampler -```{julia} -# Settings of the Hamiltonian Monte Carlo (HMC) sampler. -iterations = 1000 -ϵ = 0.05 -τ = 10 - -# Start sampling. -chain = sample(coinflip(data), HMC(ϵ, τ), iterations, progress=false) -``` - - -## Plot a summary -```{julia} -# Plot a summary of the sampling process for the parameter p, i.e. the probability of heads in a coin. -histogram(chain[:p]) -``` \ No newline at end of file From e6b4431cd439de03f935221e18bb4f96e15145c3 Mon Sep 17 00:00:00 2001 From: Penelope Yong Date: Tue, 10 Sep 2024 16:56:49 +0100 Subject: [PATCH 3/3] Remove unneeded detail from getting started --- tutorials/docs-00-getting-started/index.qmd | 120 +++----------------- 1 file changed, 18 insertions(+), 102 deletions(-) diff --git a/tutorials/docs-00-getting-started/index.qmd b/tutorials/docs-00-getting-started/index.qmd index f2d949345..e8e69a904 100644 --- a/tutorials/docs-00-getting-started/index.qmd +++ b/tutorials/docs-00-getting-started/index.qmd @@ -21,6 +21,7 @@ You will need to install Julia 1.7 or greater, which you can get from [the offic Turing is officially registered in the [Julia General package registry](https://github.com/JuliaRegistries/General), which means that you can install a stable version of Turing by running the following in the Julia REPL: ```{julia} +#| eval: false #| output: false using Pkg Pkg.add("Turing") @@ -28,8 +29,6 @@ Pkg.add("Turing") ### Example usage -Here is a simple example showing Turing in action. - First, we load the Turing and StatsPlots modules. The latter is required for visualising the results. @@ -38,8 +37,18 @@ using Turing using StatsPlots ``` -Then, we define a simple Normal model with unknown mean and variance. -Here, both `x` and `y` are observed values (since they are function parameters), whereas `m` and `s²` are the parameters to be inferred. +We then specify our model, which is a simple Gaussian model with unknown mean and variance. +In mathematical notation, the model is defined as follows: + +$$\begin{align} +s^2 &\sim \text{InverseGamma}(2, 3) \\ +m &\sim \mathcal{N}(0, \sqrt{s^2}) \\ +x, y &\sim \mathcal{N}(m, s^2) +\end{align}$$ + +This translates directly into the following Turing model. +Here, both `x` and `y` are observed values, and should therefore be passed as function parameters. +`m` and `s²` are the parameters to be inferred. ```{julia} @model function gdemo(x, y) @@ -50,115 +59,22 @@ Here, both `x` and `y` are observed values (since they are function parameters), end ``` -Suppose we have some data `x` and `y` that we want to infer the mean and variance for. +Suppose we observe `x = 1.5` and `y = 2`, and want to infer the mean and variance. We can pass these data as arguments to the `gdemo` function, and run a sampler to collect the results. -In this specific example, we collect 1000 samples using the No U-Turn Sampler (NUTS) algorithm, which is a Hamiltonian Monte Carlo method. +Here, we collect 1000 samples using the No U-Turn Sampler (NUTS) algorithm. ```{julia} -chn = sample(gdemo(1.5, 2), NUTS(), 1000, progress=false) +chain = sample(gdemo(1.5, 2), NUTS(), 1000, progress=false) ``` We can plot the results: ```{julia} -plot(chn) +plot(chain) ``` and obtain summary statistics by indexing the chain: ```{julia} -mean(chn[:m]), mean(chn[:s²]) +mean(chain[:m]), mean(chain[:s²]) ``` - - -### Verifying the results - -To verify the results of this example, we can analytically solve for the posterior distribution. -Our prior is a normal-inverse gamma distribution: - -$$\begin{align} - P(s^2) &= \Gamma^{-1}(s^2 | \alpha = 2, \beta = 3) \\ - P(m) &= \mathcal{N}(m | \mu = 0, \sigma^2 = s^2) \\ -\implies P(m, s^2) &= \text{N-}\Gamma^{-1}(m, s^2 | \mu = 0, \lambda = 1, \alpha = 2, \beta = 3). -\end{align}$$ - -This is a conjugate prior for a normal distribution with unknown mean and variance. -So, the posterior distribution given some data $\mathbf{x}$ is also a normal-inverse gamma distribution, which we can compute analytically: - -$$P(m, s^2 | \mathbf{x}) = \text{N-}\Gamma^{-1}(m, s^2 | \mu', \lambda', \alpha', \beta').$$ - -The four updated parameters are given respectively (see for example [Wikipedia](https://en.wikipedia.org/wiki/Normal-gamma_distribution#Posterior_distribution_of_the_parameters)) by - -$$\begin{align} -\mu' &= \frac{\lambda\mu + n\bar{x}}{\lambda + n}; & -\lambda' &= \lambda + n; \\[0.3cm] -\alpha' &= \alpha + \frac{n}{2}; & -\beta' &= \beta + \frac{1}{2}\left(nv + \frac{\lambda n(\bar{x} - \mu)^2}{\lambda + n}\right), -\end{align}$$ - -where $n$ is the number of data points, $\bar{x}$ the mean of the data, and $v$ the uncorrected sample variance of the data. - -```{julia} -s² = InverseGamma(2, 3) -α, β = shape(s²), scale(s²) - -# We defined m ~ Normal(0, sqrt(s²)), and the normal-inverse gamma -# distribution is defined by m ~ Normal(µ, sqrt(s²/λ)), which lets -# us identify the following -µ, λ = 0, 1 - -data = [1.5, 2] -x_bar = mean(data) -n = length(data) -v = var(data; corrected=false) - -µ´ = ((λ * µ) + (n * x_bar)) / (λ + n) -λ´ = λ + n -α´ = α + n / 2 -β´ = β + ((n * v) + ((λ * n * (x_bar - µ)^2) / (λ + n))) / 2 - -(µ´, λ´, α´, β´) -``` - -We can use these to analytically calculate, for example, the expectation values of the mean and variance parameters in the posterior distribution: $E[m] = \mu'$, and $E[s^2] = \beta'/(\alpha - - 1)$. - -```{julia} -mean_exp = µ´ -variance_exp = β´ / (α´ - 1) - -(mean_exp, variance_exp) -``` - -Alternatively, we can compare the two distributions visually (one parameter at a time). -To do so, we will directly sample from the analytical posterior distribution using the procedure described in [Wikipedia](https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution#Generating_normal-inverse-gamma_random_variates). - -```{julia} -function sample_posterior_analytic(µ´, λ´, α´, β´, iterations) - samples = [] - for i in 1:iterations - sampled_s² = rand(InverseGamma(α´, β')) - sampled_m = rand(Normal(µ´, sqrt(sampled_s² / λ´))) - samples = push!(samples, (sampled_m, sampled_s²)) - end - return samples -end - -analytical_samples = sample_posterior_analytic(µ´, λ´, α´, β´, 1000); -``` - -Comparing the distribution of the mean: - -```{julia} -density([s[1] for s in analytical_samples]; label="Posterior (Analytical)") -density!(chn[:m]; label="Posterior (HMC)") -``` - -and the distribution of the variance: - -```{julia} -density([s[2] for s in analytical_samples]; label="Posterior (Analytical)") -density!(chn[:s²]; label="Posterior (HMC)") -``` - -we see that our HMC sampler has given us a good representation of the posterior distribution.