-
-
-
## Variational Inference (ADVI)
We can also use Turing's variational inference tools to estimate the parameters of this model. See [variational inference](https://turing.ml/dev/docs/for-developers/variational_inference) for more information.
@@ -258,6 +207,7 @@ We can also use Turing's variational inference tools to estimate the parameters
```julia
using Bijectors
using Turing: Variational
+using AdvancedVI
m = bayes_nn(hcat(xs...), ts);
@@ -266,27 +216,20 @@ q = Variational.meanfield(m)
μ = randn(length(q))
ω = -1 .* ones(length(q))
-q = Variational.update(q, μ, exp.(ω));
+q = AdvancedVI.update(q, μ, exp.(ω));
-advi = ADVI(10, 1000)
+advi = ADVI(10, 5_000)
q_hat = vi(m, advi, q);
```
- ┌ Info: [ADVI] Should only be seen once: optimizer created for θ
- │ objectid(θ) = 3812708583762184342
- └ @ Turing.Variational /home/cameron/.julia/packages/Turing/cReBm/src/variational/VariationalInference.jl:204
-
-
-
```julia
samples = transpose(rand(q_hat, 5000))
-ch_vi = Chains(reshape(samples, size(samples)..., 1), ["nn_params[$$i]" for i = 1:20]);
+ch_vi = Chains(reshape(samples, size(samples)..., 1), string.(MCMCChains.namesingroup(ch, :nn_params)));
# Extract all weight and bias parameters.
-theta = ch_vi[:nn_params].value.data;
+theta = MCMCChains.group(ch_vi, :nn_params).value;
```
-
```julia
# Plot the average prediction.
plot_data()
@@ -298,13 +241,6 @@ Z = [nn_predict([x, y], theta, n_end)[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)
```
-
-
-
-
-
-
-
## Generic Bayesian Neural Networks
The below code is intended for use in more general applications, where you need to be able to change the basic network shape fluidly. The code above is highly rigid, and adapting it for other architectures would be time consuming. Currently the code below only supports networks of `Dense` layers.
@@ -366,23 +302,11 @@ end
end
end
-# Set the backend.
-Turing.setadbackend(:reverse_diff)
-
# Perform inference.
num_samples = 500
ch2 = sample(bayes_nn_general(hcat(xs...), ts, network_shape, num_params), NUTS(0.65), num_samples);
```
- ┌ Warning: `Turing.setadbackend(:reverse_diff)` is deprecated. Please use `Turing.setadbackend(:tracker)` to use `Tracker` or `Turing.setadbackend(:reversediff)` to use `ReverseDiff`. To use `ReverseDiff`, please make sure it is loaded separately with `using ReverseDiff`.
- │ caller = setadbackend(::Symbol) at ad.jl:5
- └ @ Turing.Core /home/cameron/.julia/packages/Turing/cReBm/src/core/ad.jl:5
- ┌ Info: Found initial step size
- │ ϵ = 0.2
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/cReBm/src/inference/hmc.jl:556
-
-
-
```julia
# This function makes predictions based on network shape.
function nn_predict(x, theta, num, network_shape)
@@ -390,7 +314,7 @@ function nn_predict(x, theta, num, network_shape)
end;
# Extract the θ parameters from the sampled chain.
-params2 = ch2[:θ].value.data
+params2 = MCMCChains.group(ch2, :θ).value
plot_data()
@@ -400,11 +324,4 @@ Z = [nn_predict([x, y], params2, length(ch2), network_shape)[1] for x=x_range, y
contour!(x_range, y_range, Z)
```
-
-
-
-
-
-
-
This has been an introduction to the applications of Turing and Flux in defining Bayesian neural networks.
diff --git a/tutorials/bayesian-deep-learning/Project.toml b/tutorials/bayesian-deep-learning/Project.toml
new file mode 100644
index 000000000..df324f01c
--- /dev/null
+++ b/tutorials/bayesian-deep-learning/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+AdvancedVI = "b5ca4192-6429-45e5-a2d9-87aec30a685c"
+Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
+Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
diff --git a/tutorials/differential-equations/01_bayesian-diff-eq.jmd b/tutorials/differential-equations/01_bayesian-diff-eq.jmd
index 7f01456ce..4c70407dd 100644
--- a/tutorials/differential-equations/01_bayesian-diff-eq.jmd
+++ b/tutorials/differential-equations/01_bayesian-diff-eq.jmd
@@ -25,8 +25,6 @@ $$\frac{dx}{dt} = (\alpha - \beta y)x$$
$$\frac{dy}{dt} = (\delta x - \gamma)y$$
-
-
```julia
function lotka_volterra(du,u,p,t)
x, y = u
@@ -41,30 +39,13 @@ sol = solve(prob,Tsit5())
plot(sol)
```
-
-
-
-
-
-
-
We'll generate the data to use for the parameter estimation from simulation.
With the `saveat` [argument](https://docs.sciml.ai/latest/basics/common_solver_opts/) we specify that the solution is stored only at `0.1` time units.
-
```julia
odedata = Array(solve(prob,Tsit5(),saveat=0.1))
```
-
-
-
- 2×101 Array{Float64,2}:
- 1.0 1.03981 1.05332 1.03247 0.972908 … 0.133965 0.148601 0.165247
- 1.0 1.22939 1.52387 1.88714 2.30908 0.476902 0.450153 0.426924
-
-
-
## Fitting Lotka-Volterra with DiffEqBayes
[DiffEqBayes.jl](https://github.com/SciML/DiffEqBayes.jl) is a high level package that set of extension functionality for estimating the parameters of differential equations using Bayesian methods. It allows the choice of using CmdStan.jl, Turing.jl, DynamicHMC.jl and ApproxBayes.jl to perform a Bayesian estimation of a differential equation problem specified via the DifferentialEquations.jl interface. You can read the [docs](https://docs.sciml.ai/latest/analysis/parameter_estimation/#Bayesian-Methods-1) for an understanding of the available functionality.
@@ -77,145 +58,18 @@ priors = [truncated(Normal(1.5,0.5),0.5,2.5),truncated(Normal(1.2,0.5),0,2),trun
bayesian_result_turing = turing_inference(prob,Tsit5(),t,odedata,priors,num_samples=10_000)
```
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Info: Found initial step size
- │ ϵ = 0.00625
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/GMBTf/src/inference/hmc.jl:629
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
-
-
-
-
-
- Object of type Chains, with data of type 9000×17×1 Array{Float64,3}
-
- Iterations = 1:9000
- Thinning interval = 1
- Chains = 1
- Samples per chain = 9000
- internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
- parameters = theta[1], theta[2], theta[3], theta[4], σ[1]
-
- 2-element Array{ChainDataFrame,1}
-
- Summary Statistics
- parameters mean std naive_se mcse ess r_hat
- ────────── ────── ────── ──────── ────── ───────── ──────
- theta[1] 2.3263 0.1073 0.0011 0.0021 2202.3643 1.0000
- theta[2] 1.5434 0.0957 0.0010 0.0019 2575.4033 1.0002
- theta[3] 3.1259 0.1983 0.0021 0.0031 4127.1344 1.0000
- theta[4] 1.8356 0.0827 0.0009 0.0017 2189.2825 1.0000
- σ[1] 0.8569 0.0436 0.0005 0.0005 6856.5421 0.9999
-
- Quantiles
- parameters 2.5% 25.0% 50.0% 75.0% 97.5%
- ────────── ────── ────── ────── ────── ──────
- theta[1] 2.1185 2.2428 2.3337 2.4169 2.4916
- theta[2] 1.3655 1.4750 1.5422 1.6075 1.7367
- theta[3] 2.7571 2.9893 3.1166 3.2546 3.5440
- theta[4] 1.6902 1.7708 1.8307 1.9006 1.9868
- σ[1] 0.7755 0.8266 0.8551 0.8847 0.9484
-
-
-
-
The estimated parameters are clearly very close to the desired parameter values. We can also check that the chains have converged in the plot.
-
```julia
plot(bayesian_result_turing)
```
-
-
-
-
-
-
-
## Direct Handling of Bayesian Estimation with Turing
You could want to do some sort of reduction with the differential equation's solution or use it in some other way as well. In those cases DiffEqBayes might not be useful. Turing and DifferentialEquations are completely composable and you can write of the differential equation inside a Turing `@model` and it will just work.
We can rewrite the Lotka Volterra parameter estimation problem with a Turing `@model` interface as below
-
```julia
Turing.setadbackend(:forwarddiff)
@@ -239,47 +93,6 @@ model = fitlv(odedata)
chain = sample(model, NUTS(.65),10000)
```
- ┌ Info: Found initial step size
- │ ϵ = 0.2
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/GMBTf/src/inference/hmc.jl:629
- [32mSampling: 100%|█████████████████████████████████████████| Time: 0:02:48[39m
-
-
-
-
-
- Object of type Chains, with data of type 9000×17×1 Array{Float64,3}
-
- Iterations = 1:9000
- Thinning interval = 1
- Chains = 1
- Samples per chain = 9000
- internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
- parameters = α, β, γ, δ, σ
-
- 2-element Array{ChainDataFrame,1}
-
- Summary Statistics
- parameters mean std naive_se mcse ess r_hat
- ────────── ────── ────── ──────── ────── ───────── ──────
- α 1.4999 0.0060 0.0001 0.0001 2341.1779 0.9999
- β 0.9999 0.0037 0.0000 0.0001 2440.6968 0.9999
- γ 3.0001 0.0047 0.0000 0.0001 4070.6419 1.0003
- δ 1.0001 0.0032 0.0000 0.0001 2324.4733 0.9999
- σ 0.0151 0.0011 0.0000 0.0000 4591.2728 0.9999
-
- Quantiles
- parameters 2.5% 25.0% 50.0% 75.0% 97.5%
- ────────── ────── ────── ────── ────── ──────
- α 1.4881 1.4960 1.4998 1.5038 1.5118
- β 0.9925 0.9975 0.9999 1.0024 1.0074
- γ 2.9911 2.9970 3.0000 3.0032 3.0095
- δ 0.9937 0.9979 1.0001 1.0022 1.0066
- σ 0.0131 0.0143 0.0150 0.0158 0.0173
-
-
-
-
## Scaling to Large Models: Adjoint Sensitivities
DifferentialEquations.jl's efficiency for large stiff models has been shown in multiple [benchmarks](https://github.com/SciML/DiffEqBenchmarks.jl). To learn more about how to optimize solving performance for stiff problems you can take a look at the [docs](https://docs.sciml.ai/latest/tutorials/advanced_ode_example/).
@@ -297,7 +110,9 @@ All we had to do is switch the AD backend to one of the adjoint-compatible backe
```julia
+using Zygote
Turing.setadbackend(:zygote)
+
@model function fitlv(data)
σ ~ InverseGamma(2, 3)
α ~ truncated(Normal(1.5,0.5),0.5,2.5)
@@ -315,50 +130,8 @@ model = fitlv(odedata)
chain = sample(model, NUTS(.65),1000)
```
- ┌ Info: Found initial step size
- │ ϵ = 0.2
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/GMBTf/src/inference/hmc.jl:629
- [32mSampling: 100%|█████████████████████████████████████████| Time: 0:10:42[39m
-
-
-
-
-
- Object of type Chains, with data of type 500×17×1 Array{Float64,3}
-
- Iterations = 1:500
- Thinning interval = 1
- Chains = 1
- Samples per chain = 500
- internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
- parameters = α, β, γ, δ, σ
-
- 2-element Array{ChainDataFrame,1}
-
- Summary Statistics
- parameters mean std naive_se mcse ess r_hat
- ────────── ────── ────── ──────── ────── ──────── ──────
- α 1.4997 0.0052 0.0002 0.0003 201.5277 1.0046
- β 0.9999 0.0033 0.0001 0.0001 219.1974 1.0027
- γ 3.0003 0.0047 0.0002 0.0003 290.3332 1.0014
- δ 1.0002 0.0029 0.0001 0.0002 210.0807 1.0046
- σ 0.0151 0.0010 0.0000 0.0001 246.6502 1.0017
-
- Quantiles
- parameters 2.5% 25.0% 50.0% 75.0% 97.5%
- ────────── ────── ────── ────── ────── ──────
- α 1.4892 1.4962 1.5002 1.5030 1.5108
- β 0.9934 0.9978 1.0000 1.0019 1.0066
- γ 2.9910 2.9971 3.0002 3.0039 3.0084
- δ 0.9943 0.9983 1.0000 1.0021 1.0060
- σ 0.0131 0.0143 0.0151 0.0158 0.0172
-
-
-
-
Now we can exercise control of the sensitivity analysis method that is used by using the `sensealg` keyword argument. Let's choose the `InterpolatingAdjoint` from the available AD [methods](https://docs.sciml.ai/latest/analysis/sensitivity/#Sensitivity-Algorithms-1) and enable a compiled ReverseDiff vector-Jacobian product:
-
```julia
@model function fitlv(data)
σ ~ InverseGamma(2, 3)
@@ -377,98 +150,6 @@ model = fitlv(odedata)
@time chain = sample(model, NUTS(.65),1000)
```
- ┌ Info: Found initial step size
- │ ϵ = 0.2
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/GMBTf/src/inference/hmc.jl:629
- [32mSampling: 11%|████▍ | ETA: 0:06:27[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 13%|█████▍ | ETA: 0:05:58[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 15%|██████▎ | ETA: 0:05:27[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 21%|████████▌ | ETA: 0:04:20[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 23%|█████████▎ | ETA: 0:04:03[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 24%|██████████ | ETA: 0:03:48[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 28%|███████████▌ | ETA: 0:03:27[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 29%|███████████▊ | ETA: 0:03:24[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 29%|████████████ | ETA: 0:03:20[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 36%|███████████████ | ETA: 0:02:45[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 37%|███████████████▏ | ETA: 0:02:44[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 39%|████████████████ | ETA: 0:02:36[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 46%|██████████████████▉ | ETA: 0:02:08[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 48%|███████████████████▊ | ETA: 0:02:03[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 49%|████████████████████▏ | ETA: 0:02:01[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 50%|████████████████████▎ | ETA: 0:02:00[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 100%|█████████████████████████████████████████| Time: 0:03:32[39m
-
-
- 225.663919 seconds (1.41 G allocations: 66.216 GiB, 5.25% gc time)
-
-
-
-
-
- Object of type Chains, with data of type 500×17×1 Array{Float64,3}
-
- Iterations = 1:500
- Thinning interval = 1
- Chains = 1
- Samples per chain = 500
- internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
- parameters = α, β, γ, δ, σ
-
- 2-element Array{ChainDataFrame,1}
-
- Summary Statistics
- parameters mean std naive_se mcse ess r_hat
- ────────── ────── ────── ──────── ────── ──────── ──────
- α 0.9122 0.2810 0.0126 0.0152 211.4497 0.9992
- β 1.8499 0.1141 0.0051 0.0055 302.7650 1.0018
- γ 2.5879 0.3299 0.0148 0.0228 307.5110 0.9997
- δ 0.1259 0.0221 0.0010 0.0007 219.5371 1.0006
- σ 0.8746 0.0437 0.0020 0.0017 342.6660 1.0008
-
- Quantiles
- parameters 2.5% 25.0% 50.0% 75.0% 97.5%
- ────────── ────── ────── ────── ────── ──────
- α 0.5060 0.6920 0.8932 1.0874 1.5467
- β 1.5810 1.7796 1.8709 1.9437 1.9873
- γ 1.9519 2.3707 2.5999 2.8158 3.1966
- δ 0.0843 0.1103 0.1245 0.1410 0.1704
- σ 0.7984 0.8444 0.8722 0.9044 0.9651
-
-
-
-
For more examples of adjoint usage on large parameter models, consult the [DiffEqFlux documentation](https://diffeqflux.sciml.ai/dev/)
## Including Process Noise: Estimation of Stochastic Differential Equations
@@ -483,7 +164,6 @@ $$dx = (\alpha - \beta y)xdt + \phi_1 xdW_1$$
$$dy = (\delta x - \gamma)ydt + \phi_2 ydW_2$$
-
```julia
function lotka_volterra_noise(du,u,p,t)
du[1] = p[5]*u[1]
@@ -493,18 +173,8 @@ p = [1.5, 1.0, 3.0, 1.0, 0.3, 0.3]
prob = SDEProblem(lotka_volterra,lotka_volterra_noise,u0,(0.0,10.0),p)
```
-
-
-
- [36mSDEProblem[0m with uType [36mArray{Float64,1}[0m and tType [36mFloat64[0m. In-place: [36mtrue[0m
- timespan: (0.0, 10.0)
- u0: [1.0, 1.0]
-
-
-
Solving it repeatedly confirms the randomness of the solution
-
```julia
sol = solve(prob,saveat=0.01)
p1 = plot(sol)
@@ -515,63 +185,29 @@ p3 = plot(sol)
plot(p1,p2,p3)
```
-
-
-
-
-
-
-
With the `MonteCarloSummary` it is easy to summarize the results from multiple runs through the `EnsembleProblem` interface, here we run the problem for 1000 `trajectories` and visualize the summary:
-
```julia
sol = solve(EnsembleProblem(prob),SRIW1(),saveat=0.1,trajectories=500)
summ = MonteCarloSummary(sol)
plot(summ)
```
-
-
-
-
-
-
-
Get data from the means to fit:
-
```julia
using DiffEqBase.EnsembleAnalysis
averagedata = Array(timeseries_steps_mean(sol))
```
-
-
-
- 2×101 Array{Float64,2}:
- 1.0 1.04218 1.05885 1.03187 0.967422 … 0.190811 0.197071 0.203714
- 1.0 1.22803 1.5283 1.89036 2.30967 1.16424 1.11006 1.07984
-
-
-
Now fit the means with Turing.
We will utilize multithreading with the [`EnsembleProblem`](https://docs.sciml.ai/stable/tutorials/sde_example/#Ensemble-Simulations-1) interface to speed up the SDE parameter estimation.
-
```julia
Threads.nthreads()
```
-
-
-
- 16
-
-
-
-
```julia
Turing.setadbackend(:forwarddiff)
@@ -596,63 +232,4 @@ end;
model = fitlv(averagedata)
chain = sample(model, NUTS(.65),500)
-```
-
- ┌ Info: Found initial step size
- │ ϵ = 0.2
- └ @ Turing.Inference /home/cameron/.julia/packages/Turing/GMBTf/src/inference/hmc.jl:629
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 0%|▏ | ETA: 0:03:49[39m┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- ┌ Warning: The current proposal will be rejected due to numerical error(s).
- │ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
- └ @ AdvancedHMC /home/cameron/.julia/packages/AdvancedHMC/WJCQA/src/hamiltonian.jl:47
- [32mSampling: 100%|█████████████████████████████████████████| Time: 2:33:35[39m
-
-
-
-
-
- Object of type Chains, with data of type 250×19×1 Array{Float64,3}
-
- Iterations = 1:250
- Thinning interval = 1
- Chains = 1
- Samples per chain = 250
- internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
- parameters = α, β, γ, δ, σ, ϕ1, ϕ2
-
- 2-element Array{ChainDataFrame,1}
-
- Summary Statistics
- parameters mean std naive_se mcse ess r_hat
- ────────── ────── ────── ──────── ────── ────── ──────
- α 1.6255 0.0000 0.0000 0.0000 2.0325 2.5501
- β 1.1163 0.0000 0.0000 0.0000 2.0325 Inf
- γ 3.2056 0.0000 0.0000 0.0000 2.0325 0.9960
- δ 0.9268 0.0000 0.0000 0.0000 2.0325 2.9880
- σ 0.0669 0.0000 0.0000 0.0000 2.0325 1.1011
- ϕ1 0.2329 0.0000 0.0000 0.0000 2.0325 3.2549
- ϕ2 0.2531 0.0000 0.0000 0.0000 2.0325 0.9960
-
- Quantiles
- parameters 2.5% 25.0% 50.0% 75.0% 97.5%
- ────────── ────── ────── ────── ────── ──────
- α 1.6255 1.6255 1.6255 1.6255 1.6255
- β 1.1163 1.1163 1.1163 1.1163 1.1163
- γ 3.2056 3.2056 3.2056 3.2056 3.2056
- δ 0.9268 0.9268 0.9268 0.9268 0.9268
- σ 0.0669 0.0669 0.0669 0.0669 0.0669
- ϕ1 0.2329 0.2329 0.2329 0.2329 0.2329
- ϕ2 0.2531 0.2531 0.2531 0.2531 0.2531
-
-
-
-
-
-```julia
-
```
diff --git a/tutorials/differential-equations/Project.toml b/tutorials/differential-equations/Project.toml
new file mode 100644
index 000000000..4330259fe
--- /dev/null
+++ b/tutorials/differential-equations/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DiffEqBayes = "ebbdde9d-f333-5424-9be2-dbf1e9acfb5e"
+DiffEqSensitivity = "41bf760c-e81c-5289-8e54-58b1f1f8abe2"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
diff --git a/tutorials/graphical-models/01_hidden-markov-model.jmd b/tutorials/graphical-models/01_hidden-markov-model.jmd
index 2dc5f9057..bbe0bccd6 100644
--- a/tutorials/graphical-models/01_hidden-markov-model.jmd
+++ b/tutorials/graphical-models/01_hidden-markov-model.jmd
@@ -36,13 +36,6 @@ N = length(y); K = 3;
plot(y, xlim = (0,15), ylim = (-1,5), size = (500, 250))
```
-
-
-
-
-
-
-
We can see that we have three states, one for each height of the plot (1, 2, 3). This height is also our emission parameter, so state one produces a value of one, state two produces a value of two, and so on.
Ultimately, we would like to understand three major parameters:
@@ -53,9 +46,9 @@ Ultimately, we would like to understand three major parameters:
With this in mind, let's set up our model. We are going to use some of our knowledge as modelers to provide additional information about our system. This takes the form of the prior on our emission parameter.
-\$\$
-m_i \sim Normal(i, 0.5), \space m = \{1,2,3\}
-\$\$
+$$
+m_i \sim \mathrm{Normal}(i, 0.5) \quad \text{where} \quad m = \{1,2,3\}
+$$
Simply put, this says that we expect state one to emit values in a Normally distributed manner, where the mean of each state's emissions is that state's value. The variance of 0.5 helps the model converge more quickly — consider the case where we have a variance of 1 or 2. In this case, the likelihood of observing a 2 when we are in state 1 is actually quite high, as it is within a standard deviation of the true emission value. Applying the prior that we are likely to be tightly centered around the mean prevents our model from being too confused about the state that is generating our observations.
@@ -120,8 +113,8 @@ The code below generates an animation showing the graph of the data above, and t
using StatsPlots
# Extract our m and s parameters from the chain.
-m_set = c[:m].value.data
-s_set = c[:s].value.data
+m_set = MCMCChains.group(c, :m).value
+s_set = MCMCChains.group(c, :s).value
# Iterate through the MCMC samples.
Ns = 1:length(c)
@@ -130,7 +123,7 @@ Ns = 1:length(c)
animation = @animate for i in Ns
m = m_set[i, :];
s = Int.(s_set[i,:]);
- emissions = collect(skipmissing(m[s]))
+ emissions = m[s]
p = plot(y, c = :red,
size = (500, 250),
@@ -139,60 +132,32 @@ animation = @animate for i in Ns
legend = :topright, label = "True data",
xlim = (0,15),
ylim = (-1,5));
- plot!(emissions, color = :blue, label = "Sample $$N")
-end every 10;
+ plot!(emissions, color = :blue, label = "Sample $N")
+end every 10
```
-
-
Looks like our model did a pretty good job, but we should also check to make sure our chain converges. A quick check is to examine whether the diagonal (representing the probability of remaining in the current state) of the transition matrix appears to be stationary. The code below extracts the diagonal and shows a traceplot of each persistence probability.
```julia
# Index the chain with the persistence probabilities.
-subchain = c[:,["T[$$i][$$i]" for i in 1:K],:]
+subchain = MCMCChains.group(c, :T)
+# TODO: This doesn't work anymore. Note sure what it was originally doing
# Plot the chain.
-plot(subchain,
+plot(
+ subchain,
colordim = :parameter,
seriestype=:traceplot,
title = "Persistence Probability",
legend=:right
- )
+)
```
-
-
-
-
-
-
-
A cursory examination of the traceplot above indicates that at least `T[3,3]` and possibly `T[2,2]` have converged to something resembling stationary. `T[1,1]`, on the other hand, has a slight "wobble", and seems less consistent than the others. We can use the diagnostic functions provided by [MCMCChain](https://github.com/TuringLang/MCMCChain.jl) to engage in some formal tests, like the Heidelberg and Welch diagnostic:
-
```julia
-heideldiag(c[:T])
+heideldiag(MCMCChains.group(c, :T))
```
-
-
-
- 1-element Array{ChainDataFrame{NamedTuple{(:parameters, Symbol("Burn-in"), :Stationarity, Symbol("p-value"), :Mean, :Halfwidth, :Test),Tuple{Array{String,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1}}}},1}:
- Heidelberger and Welch Diagnostic - Chain 1
- parameters Burn-in Stationarity p-value Mean Halfwidth Test
- ────────── ─────── ──────────── ─────── ────── ───────── ──────
- T[1][1] 50.0000 0.0000 0.0001 0.5329 0.0063 1.0000
- T[1][2] 50.0000 0.0000 0.0189 0.1291 0.0043 1.0000
- T[1][3] 50.0000 0.0000 0.0230 0.3381 0.0032 1.0000
- T[2][1] 30.0000 1.0000 0.2757 0.0037 0.0000 1.0000
- T[2][2] 0.0000 1.0000 0.1689 0.0707 0.0022 1.0000
- T[2][3] 0.0000 1.0000 0.1365 0.9255 0.0022 1.0000
- T[3][1] 50.0000 0.0000 0.0454 0.4177 0.0147 1.0000
- T[3][2] 40.0000 1.0000 0.0909 0.2549 0.0080 1.0000
- T[3][3] 50.0000 0.0000 0.0098 0.3274 0.0067 1.0000
-
-
-
-
The p-values on the test suggest that we cannot reject the hypothesis that the observed sequence comes from a stationary distribution, so we can be somewhat more confident that our transition matrix has converged to something reasonable.
diff --git a/tutorials/graphical-models/Project.toml b/tutorials/graphical-models/Project.toml
new file mode 100644
index 000000000..2704e9aad
--- /dev/null
+++ b/tutorials/graphical-models/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
diff --git a/tutorials/introduction/01_introduction.jmd b/tutorials/introduction/01_introduction.jmd
index 0fec4556d..452ff0a42 100644
--- a/tutorials/introduction/01_introduction.jmd
+++ b/tutorials/introduction/01_introduction.jmd
@@ -53,18 +53,6 @@ data = rand(Bernoulli(p_true), last(Ns))
data[1:5]
```
-
-
-
- 5-element Array{Bool,1}:
- 1
- 0
- 1
- 1
- 0
-
-
-
After flipping all our coins, we want to set a prior belief about what we think the distribution of coin flips look like. In this case, we are going to choose a common prior distribution called the [Beta](https://en.wikipedia.org/wiki/Beta_distribution) distribution.
@@ -81,11 +69,11 @@ For the mathematically inclined, the `Beta` distribution is updated by adding ea
This works because mean of the `Beta` distribution is defined as the following:
-\$\$ \text{E}[\text{Beta}] = \dfrac{\alpha}{\alpha+\beta} \$\$
+$$\text{E}[\text{Beta}] = \dfrac{\alpha}{\alpha+\beta}$$
Which is 0.5 when $$\alpha = \beta$$, as we expect for a large enough number of coin flips. As we increase the number of samples, our variance will also decrease, such that the distribution will reflect less uncertainty about the probability of receiving a heads. The definition of the variance for the `Beta` distribution is the following:
-\$\$ \text{var}[\text{Beta}] = \dfrac{\alpha\beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)} \$\$
+$$\text{var}[\text{Beta}] = \dfrac{\alpha\beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}$$
The intuition about this definition is that the variance of the distribution will approach 0 with more and more samples, as the denominator will grow faster than will the numerator. More samples means less variance.
@@ -107,23 +95,16 @@ animation = @gif for (i, N) in enumerate(Ns)
# Plotting
plot(updated_belief,
size = (500, 250),
- title = "Updated belief after $$N observations",
+ title = "Updated belief after $N observations",
xlabel = "probability of heads",
ylabel = "",
legend = nothing,
xlim = (0,1),
fill=0, α=0.3, w=3)
vline!([p_true])
-end;
+end
```
- ┌ Info: Saved animation to
- │ fn = /home/cameron/code/TuringTutorials/tmp.gif
- └ @ Plots /home/cameron/.julia/packages/Plots/Xnzc7/src/animation.jl:104
-
-
-
-
The animation above shows that with increasing evidence our belief about the probability of heads in a coin flip slowly adjusts towards the true value. The orange line in the animation represents the true probability of seeing heads on a single coin flip, while the mode of the distribution shows what the model believes the probability of a heads is given the evidence it has seen.
### Coin Flipping With Turing
@@ -151,7 +132,7 @@ First, we define the coin-flip model using Turing.
@model coinflip(y) = begin
# Our prior belief about the probability of heads in a coin.
- p ~ Beta(1, 1)
+ p ~ Beta(1, 1)
# The number of observations.
N = length(y)
@@ -184,13 +165,6 @@ p_summary = chain[:p]
plot(p_summary, seriestype = :histogram)
```
-
-
-
-
-
-
-
Now we can build our plot:
@@ -212,11 +186,4 @@ plot!(p, range(0, stop = 1, length = 100), pdf.(Ref(updated_belief), range(0, st
vline!(p, [p_true], label = "True probability", c = :red)
```
-
-
-
-
-
-
-
As we can see, the Turing model closely approximates the true probability. Hopefully this tutorial has provided an easy-to-follow, yet informative introduction to Turing's simpler applications. More advanced usage will be demonstrated in later tutorials.
diff --git a/tutorials/introduction/02_gaussian-mixture-model.jmd b/tutorials/introduction/02_gaussian-mixture-model.jmd
index 4ed0c4110..e444183ba 100644
--- a/tutorials/introduction/02_gaussian-mixture-model.jmd
+++ b/tutorials/introduction/02_gaussian-mixture-model.jmd
@@ -27,31 +27,24 @@ x = mapreduce(c -> rand(MvNormal([μs[c], μs[c]], 1.), N), hcat, 1:2)
scatter(x[1,:], x[2,:], legend = false, title = "Synthetic Dataset")
```
-
-
-
-
-
-
-
## Gaussian Mixture Model in Turing
To cluster the data points shown above, we use a model that consists of two mixture components (clusters) and assigns each datum to one of the components. The assignment thereof determines the distribution that the data point is generated from.
In particular, in a Bayesian Gaussian mixture model with $$1 \leq k \leq K$$ components for 1-D data each data point $$x_i$$ with $$1 \leq i \leq N$$ is generated according to the following generative process.
First we draw the parameters for each cluster, i.e. in our example we draw location of the distributions from a Normal:
-\$\$
-\mu_k \sim Normal() \, , \; \forall k \\
-\$\$
+$$
+\mu_k \sim \mathrm{Normal}() \, , \; \forall k
+$$
and then draw mixing weight for the $$K$$ clusters from a Dirichlet distribution, i.e.
-\$\$
- w \sim Dirichlet(K, \alpha) \, . \\
-\$\$
+$$
+ w \sim \mathrm{Dirichlet}(K, \alpha) \, .
+$$
After having constructed all the necessary model parameters, we can generate an observation by first selecting one of the clusters and then drawing the datum accordingly, i.e.
-\$\$
- z_i \sim Categorical(w) \, , \; \forall i \\
- x_i \sim Normal(\mu_{z_i}, 1.) \, , \; \forall i
-\$\$
+$$
+ z_i \sim \mathrm{Categorical}(w) \, , \; \forall i \\
+ x_i \sim \mathrm{Normal}(\mu_{z_i}, 1.) \, , \; \forall i
+$$
For more details on Gaussian mixture models, we refer to Christopher M. Bishop, *Pattern Recognition and Machine Learning*, Section 9.
@@ -60,21 +53,9 @@ For more details on Gaussian mixture models, we refer to Christopher M. Bishop,
using Turing, MCMCChains
# Turn off the progress monitor.
-Turing.turnprogress(false)
+Turing.turnprogress(false);
```
- ┌ Info: [Turing]: progress logging is disabled globally
- └ @ Turing /home/cameron/.julia/packages/Turing/cReBm/src/Turing.jl:22
-
-
-
-
-
- false
-
-
-
-
```julia
@model GaussianMixtureModel(x) = begin
@@ -107,13 +88,6 @@ Turing.turnprogress(false)
end
```
-
-
-
- ##GaussianMixtureModel#361 (generic function with 2 methods)
-
-
-
After having specified the model in Turing, we can construct the model function and run a MCMC simulation to obtain assignments of the data points.
@@ -139,17 +113,10 @@ In particular, in this example we consider the sample values of the location par
```julia
-ids = findall(map(name -> occursin("μ", name), names(tchain)));
+ids = findall(map(name -> occursin("μ", string(name)), names(tchain)));
p=plot(tchain[:, ids, :], legend=true, labels = ["Mu 1" "Mu 2"], colordim=:parameter)
```
-
-
-
-
-
-
-
You'll note here that it appears the location means are switching between chains. We will address this in future tutorials. For those who are keenly interested, see [this](https://mc-stan.org/users/documentation/case-studies/identifying_mixture_models.html) article on potential solutions.
For the moment, we will just use the first chain to ensure the validity of our inference.
@@ -173,44 +140,23 @@ function predict(x, y, w, μ)
end
```
-
-
-
- predict (generic function with 1 method)
-
-
-
-
```julia
contour(range(-5, stop = 3), range(-6, stop = 2),
- (x, y) -> predict(x, y, [0.5, 0.5], [mean(tchain[:μ1].value), mean(tchain[:μ2].value)])
+ (x, y) -> predict(x, y, [0.5, 0.5], [mean(tchain[:μ1]), mean(tchain[:μ2])])
)
scatter!(x[1,:], x[2,:], legend = false, title = "Synthetic Dataset")
```
-
-
-
-
-
-
-
## Infered Assignments
Finally, we can inspect the assignments of the data points infered using Turing. As we can see, the dataset is partitioned into two distinct groups.
```julia
-assignments = collect(skipmissing(mean(tchain[:k].value, dims=1).data))
+# TODO: is there a better way than this icky `.nt.mean` stuff?
+assignments = mean(MCMCChains.group(tchain, :k)).nt.mean
scatter(x[1,:], x[2,:],
legend = false,
title = "Assignments on Synthetic Dataset",
zcolor = assignments)
```
-
-
-
-
-
-
-
diff --git a/tutorials/introduction/Project.toml b/tutorials/introduction/Project.toml
new file mode 100644
index 000000000..c3d3e22a9
--- /dev/null
+++ b/tutorials/introduction/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
diff --git a/tutorials/non-parameteric/01_infinite-mixture-model.jmd b/tutorials/non-parameteric/01_infinite-mixture-model.jmd
index 310cad70a..5ac4ef9d8 100644
--- a/tutorials/non-parameteric/01_infinite-mixture-model.jmd
+++ b/tutorials/non-parameteric/01_infinite-mixture-model.jmd
@@ -23,12 +23,12 @@ The generative process of such a model can be written as:
$$
\begin{align}
-\pi_1 &\sim Beta(a, b) \\
+\pi_1 &\sim \mathrm{Beta}(a, b) \\
\pi_2 &= 1-\pi_1 \\
-\mu_1 &\sim Normal(\mu_0, \Sigma_0) \\
-\mu_2 &\sim Normal(\mu_0, \Sigma_0) \\
-z_i &\sim Categorical(\pi_1, \pi_2) \\
-x_i &\sim Normal(\mu_{z_i}, \Sigma)
+\mu_1 &\sim \mathrm{Normal}(\mu_0, \Sigma_0) \\
+\mu_2 &\sim \mathrm{Normal}(\mu_0, \Sigma_0) \\
+z_i &\sim \mathrm{Categorical}(\pi_1, \pi_2) \\
+x_i &\sim \mathrm{Normal}(\mu_{z_i}, \Sigma)
\end{align}
$$
@@ -36,41 +36,32 @@ where $$\pi_1, \pi_2$$ are the mixing weights of the mixture model, i.e. $$\pi_1
We can implement this model in Turing for 1D data as follows:
-
```julia
-@model two_model(x) = begin
-
+@model function two_model(x)
# Hyper-parameters
μ0 = 0.0
σ0 = 1.0
# Draw weights.
- π1 ~ Beta(1,1)
+ π1 ~ Beta(1,1)
π2 = 1-π1
# Draw locations of the components.
- μ1 ~ Normal(μ0, σ0)
- μ2 ~ Normal(μ0, σ0)
+ μ1 ~ Normal(μ0, σ0)
+ μ2 ~ Normal(μ0, σ0)
# Draw latent assignment.
- z ~ Categorical([π1, π2])
+ z ~ Categorical([π1, π2])
# Draw observation from selected component.
if z == 1
- x ~ Normal(μ1, 1.0)
+ x ~ Normal(μ1, 1.0)
else
- x ~ Normal(μ2, 1.0)
+ x ~ Normal(μ2, 1.0)
end
end
```
-
-
-
- DynamicPPL.ModelGen{var"###generator#282",(:x,),(),Tuple{}}(##generator#282, NamedTuple())
-
-
-
#### Finite Mixture Model
If we have more than two components, this model can elegantly be extend using a Dirichlet distribution as prior for the mixing weights $$\pi_1, \dots, \pi_K$$. Note that the Dirichlet distribution is the multivariate generalization of the beta distribution. The resulting model can be written as:
@@ -78,9 +69,9 @@ If we have more than two components, this model can elegantly be extend using a
$$
\begin{align}
(\pi_1, \dots, \pi_K) &\sim Dirichlet(K, \alpha) \\
-\mu_k &\sim Normal(\mu_0, \Sigma_0), \;\; \forall k \\
+\mu_k &\sim \mathrm{Normal}(\mu_0, \Sigma_0), \;\; \forall k \\
z &\sim Categorical(\pi_1, \dots, \pi_K) \\
-x &\sim Normal(\mu_z, \Sigma)
+x &\sim \mathrm{Normal}(\mu_z, \Sigma)
\end{align}
$$
@@ -101,9 +92,9 @@ We now will utilize the fact that one can integrate out the mixing weights in a
In fact, if the mixing weights are integrated out, the conditional prior for the latent variable $$z$$ is given by:
-\$\$
-p(z_i = k \mid z_{\not i}, \alpha) = \frac{n_k + \alpha/K}{N - 1 + \alpha}
-\$\$
+$$
+p(z_i = k \mid z_{\not i}, \alpha) = \frac{n_k + \alpha K}{N - 1 + \alpha}
+$$
where $$z_{\not i}$$ are the latent assignments of all observations except observation $$i$$. Note that we use $$n_k$$ to denote the number of observations at component $$k$$ excluding observation $$i$$. The parameter $$\alpha$$ is the concentration parameter of the Dirichlet distribution used as prior over the mixing weights.
@@ -113,15 +104,15 @@ To obtain the Chinese restaurant process construction, we can now derive the con
For $$n_k > 0$$ we obtain:
-\$\$
-p(z_i = k \mid z_{\not i}, \alpha) = \frac{n_k}{N - 1 + \alpha}
-\$\$
+$$
+p(z_i = k \mid z_{\not i}, \alpha) = \frac{n_k}{N - 1 + \alpha}
+$$
and for all infinitely many clusters that are empty (combined) we get:
-\$\$
+$$
p(z_i = k \mid z_{\not i}, \alpha) = \frac{\alpha}{N - 1 + \alpha}
-\$\$
+$$
Those equations show that the conditional prior for component assignments is proportional to the number of such observations, meaning that the Chinese restaurant process has a rich get richer property.
@@ -159,21 +150,14 @@ using Plots
# Plot the cluster assignments over time
@gif for i in 1:Nmax
scatter(collect(1:i), z[1:i], markersize = 2, xlabel = "observation (i)", ylabel = "cluster (k)", legend = false)
-end;
+end
```
- ┌ Info: Saved animation to
- │ fn = /home/cameron/code/TuringTutorials/tmp.gif
- └ @ Plots /home/cameron/.julia/packages/Plots/Xnzc7/src/animation.jl:104
-
-
-
-
Further, we can see that the number of clusters is logarithmic in the number of observations and data points. This is a side-effect of the "rich-get-richer" phenomenon, i.e. we expect large clusters and thus the number of clusters has to be smaller than the number of observations.
-\$\$
-E[K \mid N] \approx \alpha \cdot log \big(1 + \frac{N}{\alpha}\big)
-\$\$
+$$
+\mathbb{E}[K \mid N] \approx \alpha \cdot log \big(1 + \frac{N}{\alpha}\big)
+$$
We can see from the equation that the concentration parameter $$\alpha$$ allows us to control the number of clusters formed *a priori*.
@@ -181,8 +165,7 @@ In Turing we can implement an infinite Gaussian mixture model using the Chinese
```julia
-@model infiniteGMM(x) = begin
-
+@model function infiniteGMM(x)
# Hyper-parameters, i.e. concentration parameter and parameters of H.
α = 1.0
μ0 = 0.0
@@ -207,14 +190,14 @@ In Turing we can implement an infinite Gaussian mixture model using the Chinese
nk = Vector{Int}(map(k -> sum(z .== k), 1:K))
# Draw the latent assignment.
- z[i] ~ ChineseRestaurantProcess(rpm, nk)
+ z[i] ~ ChineseRestaurantProcess(rpm, nk)
# Create a new cluster?
if z[i] > K
push!(μ, 0.0)
# Draw location of new cluster.
- μ[z[i]] ~ H
+ μ[z[i]] ~ H
end
# Draw observation.
@@ -223,13 +206,6 @@ In Turing we can implement an infinite Gaussian mixture model using the Chinese
end
```
-
-
-
- DynamicPPL.ModelGen{var"###generator#800",(:x,),(),Tuple{}}(##generator#800, NamedTuple())
-
-
-
We can now use Turing to infer the assignments of some data points. First, we will create some random data that comes from three clusters, with means of 0, -5, and 10.
@@ -254,60 +230,45 @@ model_fun = infiniteGMM(data);
chain = sample(model_fun, SMC(), iterations);
```
- [32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:00[39m
-
-
Finally, we can plot the number of clusters in each sample.
```julia
# Extract the number of clusters for each sample of the Markov chain.
-k = map(t -> length(unique(chain[:z].value[t,:,:])), 1:iterations);
+k = map(
+ t -> length(unique(vec(chain[t, MCMCChains.namesingroup(chain, :z), :].value))),
+ 1:iterations
+);
# Visualize the number of clusters.
plot(k, xlabel = "Iteration", ylabel = "Number of clusters", label = "Chain 1")
```
-
-
-
-
-
-
-
If we visualize the histogram of the number of clusters sampled from our posterior, we observe that the model seems to prefer 3 clusters, which is the true number of clusters. Note that the number of clusters in a Dirichlet process mixture model is not limited a priori and will grow to infinity with probability one. However, if conditioned on data the posterior will concentrate on a finite number of clusters enforcing the resulting model to have a finite amount of clusters. It is, however, not given that the posterior of a Dirichlet process Gaussian mixture model converges to the true number of clusters, given that data comes from a finite mixture model. See Jeffrey Miller and Matthew Harrison: [A simple example of Dirichlet process mixture inconsitency for the number of components](https://arxiv.org/pdf/1301.2708.pdf) for details.
-
```julia
histogram(k, xlabel = "Number of clusters", legend = false)
```
-
-
-
-
-
-
-
One issue with the Chinese restaurant process construction is that the number of latent parameters we need to sample scales with the number of observations. It may be desirable to use alternative constructions in certain cases. Alternative methods of constructing a Dirichlet process can be employed via the following representations:
Size-Biased Sampling Process
-\$\$
-j_k \sim Beta(1, \alpha) * surplus
-\$\$
+$$
+j_k \sim \mathrm{Beta}(1, \alpha) \cdot \mathrm{surplus}
+$$
Stick-Breaking Process
-\$\$
-v_k \sim Beta(1, \alpha)
-\$\$
+$$
+v_k \sim \mathrm{Beta}(1, \alpha)
+$$
Chinese Restaurant Process
-\$\$
+$$
p(z_n = k | z_{1:n-1}) \propto \begin{cases}
\frac{m_k}{n-1+\alpha}, \text{ if } m_k > 0\\\
\frac{\alpha}{n-1+\alpha}
\end{cases}
-\$\$
+$$
For more details see [this article](https://www.stats.ox.ac.uk/~teh/research/npbayes/Teh2010a.pdf).
diff --git a/tutorials/non-parameteric/Project.toml b/tutorials/non-parameteric/Project.toml
new file mode 100644
index 000000000..dc3ff2545
--- /dev/null
+++ b/tutorials/non-parameteric/Project.toml
@@ -0,0 +1,4 @@
+[deps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
diff --git a/tutorials/regression/01_logistic-regression.jmd b/tutorials/regression/01_logistic-regression.jmd
index 037ec08b2..7e4544379 100644
--- a/tutorials/regression/01_logistic-regression.jmd
+++ b/tutorials/regression/01_logistic-regression.jmd
@@ -33,17 +33,6 @@ Random.seed!(0);
Turing.turnprogress(false)
```
- ┌ Info: [Turing]: progress logging is disabled globally
- └ @ Turing /home/cameron/.julia/packages/Turing/cReBm/src/Turing.jl:22
-
-
-
-
-
- false
-
-
-
## Data Cleaning & Set Up
Now we're going to import our dataset. The first six rows of the dataset are shown below so you can get a good feel for what kind of data we have.
@@ -57,13 +46,6 @@ data = RDatasets.dataset("ISLR", "Default");
first(data, 6)
```
-
-
-
-| Default | Student | Balance | Income | |
|---|---|---|---|---|
| Categorical… | Categorical… | Float64 | Float64 | |
| 1 | No | No | 729.526 | 44361.6 |
| 2 | No | Yes | 817.18 | 12106.1 |
| 3 | No | No | 1073.55 | 31767.1 |
| 4 | No | No | 529.251 | 35704.5 |
| 5 | No | No | 785.656 | 38463.5 |
| 6 | No | Yes | 919.589 | 7491.56 |
| Balance | Income | DefaultNum | StudentNum | |
|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | |
| 1 | 729.526 | 44361.6 | 0.0 | 0.0 |
| 2 | 817.18 | 12106.1 | 0.0 | 1.0 |
| 3 | 1073.55 | 31767.1 | 0.0 | 0.0 |
| 4 | 529.251 | 35704.5 | 0.0 | 0.0 |
| 5 | 785.656 | 38463.5 | 0.0 | 0.0 |
| 6 | 919.589 | 7491.56 | 0.0 | 1.0 |
| Model | MPG | Cyl | Disp | HP | DRat | WT | QSec | VS | |
|---|---|---|---|---|---|---|---|---|---|
| String | Float64 | Int64 | Float64 | Int64 | Float64 | Float64 | Float64 | Int64 | |
| 1 | Mazda RX4 | 21.0 | 6 | 160.0 | 110 | 3.9 | 2.62 | 16.46 | 0 |
| 2 | Mazda RX4 Wag | 21.0 | 6 | 160.0 | 110 | 3.9 | 2.875 | 17.02 | 0 |
| 3 | Datsun 710 | 22.8 | 4 | 108.0 | 93 | 3.85 | 2.32 | 18.61 | 1 |
| 4 | Hornet 4 Drive | 21.4 | 6 | 258.0 | 110 | 3.08 | 3.215 | 19.44 | 1 |
| 5 | Hornet Sportabout | 18.7 | 8 | 360.0 | 175 | 3.15 | 3.44 | 17.02 | 0 |
| 6 | Valiant | 18.1 | 6 | 225.0 | 105 | 2.76 | 3.46 | 20.22 | 1 |
| MPG | Bayes | OLS | |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 19.2 | 18.3766 | 18.1265 |
| 2 | 15.0 | 6.4176 | 6.37891 |
| 3 | 16.4 | 13.9125 | 13.883 |
| 4 | 14.3 | 11.8393 | 11.7337 |
| 5 | 21.4 | 25.3622 | 25.1916 |
| 6 | 18.1 | 20.7687 | 20.672 |
| 7 | 19.7 | 16.03 | 15.8408 |
| 8 | 15.2 | 18.2903 | 18.3391 |
| 9 | 26.0 | 28.5191 | 28.4865 |
| 10 | 17.3 | 14.498 | 14.534 |
| nsneeze | alcohol_taken | nomeds_taken | product_alcohol_meds | |
|---|---|---|---|---|
| Int64 | Float64 | Float64 | Float64 | |
| 1 | 5 | 0.0 | 0.0 | 0.0 |
| 2 | 5 | 1.0 | 0.0 | 0.0 |
| 3 | 8 | 1.0 | 0.0 | 0.0 |
| 4 | 1 | 0.0 | 0.0 | 0.0 |
| 5 | 38 | 1.0 | 1.0 | 1.0 |
| SepalLength | SepalWidth | PetalLength | PetalWidth | Species | |
|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Cat… | |
| 1 | 5.6 | 2.9 | 3.6 | 1.3 | versicolor |
| 2 | 5.8 | 2.7 | 3.9 | 1.2 | versicolor |
| 3 | 5.5 | 2.3 | 4.0 | 1.3 | versicolor |
| 4 | 6.7 | 3.3 | 5.7 | 2.5 | virginica |
| 5 | 5.1 | 3.5 | 1.4 | 0.2 | setosa |
| 6 | 5.1 | 3.8 | 1.5 | 0.3 | setosa |
| 7 | 4.8 | 3.4 | 1.9 | 0.2 | setosa |
| 8 | 6.0 | 2.9 | 4.5 | 1.5 | versicolor |
| 9 | 6.9 | 3.1 | 5.4 | 2.1 | virginica |
| 10 | 5.4 | 3.9 | 1.7 | 0.4 | setosa |
| 11 | 5.0 | 3.6 | 1.4 | 0.2 | setosa |
| 12 | 5.7 | 3.0 | 4.2 | 1.2 | versicolor |
| 13 | 5.0 | 3.3 | 1.4 | 0.2 | setosa |
| 14 | 7.7 | 3.0 | 6.1 | 2.3 | virginica |
| 15 | 5.8 | 2.8 | 5.1 | 2.4 | virginica |
| 16 | 4.4 | 3.0 | 1.3 | 0.2 | setosa |
| 17 | 6.3 | 3.3 | 4.7 | 1.6 | versicolor |
| 18 | 6.0 | 2.7 | 5.1 | 1.6 | versicolor |
| 19 | 4.6 | 3.4 | 1.4 | 0.3 | setosa |
| 20 | 6.0 | 2.2 | 4.0 | 1.0 | versicolor |
| Species | Species_index | |
|---|---|---|
| Cat… | Union… | |
| 1 | versicolor | 2 |
| 2 | setosa | 1 |
| 3 | setosa | 1 |
| 4 | versicolor | 2 |
| 5 | setosa | 1 |
| 6 | virginica | 3 |
| 7 | versicolor | 2 |
| 8 | versicolor | 2 |
| 9 | setosa | 1 |
| 10 | virginica | 3 |
| 11 | setosa | 1 |
| 12 | setosa | 1 |
| 13 | versicolor | 2 |
| 14 | versicolor | 2 |
| 15 | virginica | 3 |
| 16 | versicolor | 2 |
| 17 | setosa | 1 |
| 18 | setosa | 1 |
| 19 | virginica | 3 |
| 20 | setosa | 1 |
| Model | MPG | Cyl | Disp | HP | DRat | WT | QSec | VS | |
|---|---|---|---|---|---|---|---|---|---|
| String | Float64 | Int64 | Float64 | Int64 | Float64 | Float64 | Float64 | Int64 | |
| 1 | Mazda RX4 | 21.0 | 6 | 160.0 | 110 | 3.9 | 2.62 | 16.46 | 0 |
| 2 | Mazda RX4 Wag | 21.0 | 6 | 160.0 | 110 | 3.9 | 2.875 | 17.02 | 0 |
| 3 | Datsun 710 | 22.8 | 4 | 108.0 | 93 | 3.85 | 2.32 | 18.61 | 1 |
| 4 | Hornet 4 Drive | 21.4 | 6 | 258.0 | 110 | 3.08 | 3.215 | 19.44 | 1 |
| 5 | Hornet Sportabout | 18.7 | 8 | 360.0 | 175 | 3.15 | 3.44 | 17.02 | 0 |
| 6 | Valiant | 18.1 | 6 | 225.0 | 105 | 2.76 | 3.46 | 20.22 | 1 |
| MPG | Cyl | Disp | HP | DRat | WT | QSec | VS | |
|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 116.195 | 1.01488 | 0.591245 | 0.0483133 | -0.835198 | 0.222544 | -0.307089 | -0.868028 |
| 2 | 114.295 | 1.01488 | 0.962396 | 1.4339 | 0.249566 | 0.636461 | -1.36476 | -0.868028 |
| 3 | 120.195 | 1.01488 | 1.36582 | 0.412942 | -0.966118 | 0.641571 | -0.446992 | -0.868028 |
| 4 | 128.295 | -1.22486 | -1.22417 | -1.17684 | 0.904164 | -1.31048 | 0.588295 | 1.11604 |
| 5 | 126.995 | -1.22486 | -0.890939 | -0.812211 | 1.55876 | -1.10097 | -0.642858 | -0.868028 |
| 6 | 131.395 | -1.22486 | -1.09427 | -0.491337 | 0.324377 | -1.74177 | -0.530935 | 1.11604 |