Replication of PyMC3's simple introduction example in Julia, for the purpose of my understanding

bioassay_metropolis.jl

### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 8ce57352-cf5b-11eb-2624-976564fd9572
using Plots

# ╔═╡ 22e5afa5-03a5-4c6c-9744-7417f1176839
using Distributions

# ╔═╡ f5bba1ce-8a17-4177-9d79-f5572a0e1a78
using Random

# ╔═╡ afe35a50-90ba-46e7-b9c6-43afc93639ca
md"""## Bioassay experiment

This notebook very closely mimics (replicates) the example at https://docs.pymc.io/about.html#usage-overview, with reproduced results for the inferred values of the distributions.

The difference is that this notebook uses purely the Metropolis-Hastings algorithm (naive, no adaptiveness or similar improvements).
"""

# ╔═╡ efecda10-2dde-4133-ae4e-0189fc162a59
gr(html_output_format="png")

# ╔═╡ 53ad7217-b5cc-44b0-86b0-316dcad24a22
mutable struct Model
	observed::Vector{Float64} # observed data, for calculating the likelihood
	params::Tuple # probabilistic parameters for prior distributions etc.
	likelihood::Function # likelihood function, depends on params.
end

# ╔═╡ 257a83bb-5eea-4acc-be29-dea16081e963
begin
	# data just as in the original example
	n = fill(5, 4)
	y = [0,1,3,5]
	dose = [-.86,-.3,-.05,.73]
end

# ╔═╡ e1d2f70c-3d49-492e-bd63-a9f0118ef11c
function likelihood(params)
	((amu, asig), (bmu, bsig)) = params
	alpha = rand(Normal(amu, asig))
	beta = rand(Normal(bmu, bsig))
	theta = 1 ./(1 .+ exp.( - (alpha .+ beta .* dose))) # invlogit
	loglik = sum(log.(pdf.(Binomial.(n, theta), y)))

	(loglik, (alpha, beta))
end

# ╔═╡ 2ed0bee7-abfb-4fc6-934f-b682948dd343
# to be fair: we have a lower sigma on alpha, because the naive algorithm isn't
# good enough to deal with a very wide distribution.
mod = Model(y, ((0, 1), (0, 1)), likelihood)

# ╔═╡ 9f27fe22-5970-4aae-a14b-7b7c04d795fb
function sample_model(model::Model, n=5000, tune=div(n, 5))
	(oldlik, initial) = model.likelihood(model.params)
	samples = fill(initial, n)
	uniform = Uniform(0, 1)
	
	# This is the core of Metropolis-Hastings
	for i = 1:n
		lik, new = model.likelihood(model.params)
		ratio = exp(lik-oldlik)
		unisample = rand(uniform)
		if unisample < ratio
			oldlik = lik
			samples[i] = new
			initial = new
		else
			samples[i] = initial
		end
	end
	samples[tune:end]
end

# ╔═╡ 52c971ce-1ac9-4c5a-9db9-31592195739b
# run our "study"
samples = sample_model(mod, 40000);

# ╔═╡ 19a0c0aa-3909-4696-8807-9a3314944068
# extract distributions for alpha/beta
as, bs = [e[1] for e in samples], [e[2] for e in samples];

# ╔═╡ 432c293b-367e-41ae-9cd1-7ead81457b04
begin
	histogram(as, label="alpha")
	histogram!(bs, label="beta")
end

# ╔═╡ 5867d76b-cfff-4e3e-aa00-b97575ac4729
mean(as), mean(bs)

# ╔═╡ 2e46d17c-3ddc-4662-b2e2-d7e322f8100f
std(as), std(bs)

# ╔═╡ 9880646b-58b1-47cc-96de-c8f99e967e77
plot([as bs], label=["alpha" "beta"])

# ╔═╡ Cell order:
# ╠═afe35a50-90ba-46e7-b9c6-43afc93639ca
# ╠═8ce57352-cf5b-11eb-2624-976564fd9572
# ╠═22e5afa5-03a5-4c6c-9744-7417f1176839
# ╠═f5bba1ce-8a17-4177-9d79-f5572a0e1a78
# ╠═efecda10-2dde-4133-ae4e-0189fc162a59
# ╠═53ad7217-b5cc-44b0-86b0-316dcad24a22
# ╠═257a83bb-5eea-4acc-be29-dea16081e963
# ╠═e1d2f70c-3d49-492e-bd63-a9f0118ef11c
# ╠═2ed0bee7-abfb-4fc6-934f-b682948dd343
# ╠═9f27fe22-5970-4aae-a14b-7b7c04d795fb
# ╠═52c971ce-1ac9-4c5a-9db9-31592195739b
# ╠═19a0c0aa-3909-4696-8807-9a3314944068
# ╠═432c293b-367e-41ae-9cd1-7ead81457b04
# ╠═5867d76b-cfff-4e3e-aa00-b97575ac4729
# ╠═2e46d17c-3ddc-4662-b2e2-d7e322f8100f
# ╠═9880646b-58b1-47cc-96de-c8f99e967e77