question about quap function about sr2turingpluto.jl HOT 3 CLOSED

md" ### Quadratic approximation."
map_estimate = optimize(m4_1t, MAP())
begin
	â, b̂, σ̂ = map_estimate.values
	(â = â, b̂ = b̂, σ̂ = σ̂)
end
let
	x = 30:0.1:65
	plot(x, â .+ b̂ .* x)
	scatter!(howell1.weight, howell1.height, leg=false)
end

### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 2d3d462e-d6bf-4ebd-ab34-db38f4261a11
using Pkg, DrWatson

# ╔═╡ 18b5c138-9324-4c5a-86ca-be1e825ed11e
begin
	using PlutoUI
	using Random
	using Distributions
	using StatsPlots
	using StatsBase
	using LaTeXStrings
	using CSV
	using DataFrames
	using PrettyTables
	using StatisticalRethinking
	using StatisticalRethinkingPlots
	using LinearAlgebra
	using Logging
	using AxisKeys
	using NamedDims
	using Random
	using Turing
	using ParetoSmooth
end

# ╔═╡ 58b9da10-8cbd-42e4-bddf-ec4995ac756c
md" ## Turing guide"

# ╔═╡ d3c6c774-a158-485c-acd1-e366c24df4bc
md" ##### Widen notebook."

# ╔═╡ 0f75e85c-4506-483f-b0f8-a31ee698bf04
html"""
<style>
	main {
		margin: 0 auto;
		max-width: 2000px;
    	padding-left: max(160px, 10%);
    	padding-right: max(160px, 10%);
	}
</style>
"""

# ╔═╡ cb180e21-3778-40f1-93d8-af442d1be559
md" ##### Commonly used packages in SR2TuringPluto"

# ╔═╡ a2d950f8-d780-4abe-af37-b189bd4a5e88
md"##### Setting default attributes for logging and plots."

# ╔═╡ bbb4a889-77e2-449b-9258-883e7434e641
begin
	Logging.disable_logging(Logging.Warn);
end;

# ╔═╡ 496b257c-ea62-40b2-b1c9-e92fdda1d688
begin
	howell = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	howell1 = howell[howell.age .>= 18, :]
	PRECIS(howell1)
end

# ╔═╡ 96760a8e-8561-4576-a7c7-70596fa20d33
md" ### Statistical Rethinking (SR2) models."

# ╔═╡ b83bc4d0-8c65-40a8-9187-70ef2a0ff38d
md"In the book and associated R package `rethinking`, statistical models are defined as illustrated below:

# ╔═╡ b583eb23-4045-420c-80d5-4f339ebb6ad4
md" ##### An equivalent Turing model:"

# ╔═╡ d5d66e16-487e-4214-b702-728dc3db32cc
@model function m4_1(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* weights
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

# ╔═╡ b14217e8-1f70-44cc-8d3c-528ab55e27ef
md"
!!! note

The sequence of the statements matter in Turing models!"

# ╔═╡ b3adef33-c111-4d20-b277-a5346eae9f23
begin
	m4_1t = m4_1(howell1.weight, howell1.height)
	post4_1t = sample(m4_1t, NUTS(), MCMCThreads(), 1000, 4)
	post4_1t_df = DataFrame(post4_1t)
	PRECIS(post4_1t_df)
end

# ╔═╡ 6869d583-1590-4210-b696-b63fd7857e09
plot(post4_1t)

# ╔═╡ 1f4e22a4-5fba-4331-9bf9-dd886ce45e74
begin
	x = 30.0:0.1:64.0
	preds = mean(post4_1t[:a]) .+ mean(post4_1t[:b]) .* x
	plot(x, preds, color=:darkblue, label="Regression line")
	scatter!(howell1.weight, howell1.height, marker=:cross, markersize=3,
		color=:darkred, label="Observations")
end

# ╔═╡ cbe1028d-829d-4ed4-840d-d74f151f55bb
CHNS(post4_1t[[:a, :b, :σ, :lp, :log_density]])

# ╔═╡ 75ee3adf-2988-4f42-a0b5-21abe5829120
md" ### Priors of Turing models."

# ╔═╡ 63db6c9f-ba40-41a2-975d-798d9f2cbeab
let
	priors4_1t = sample(m4_1t, Prior(), 1000)
	PRECIS(DataFrame(priors4_1t))
end

# ╔═╡ 534edcc4-2dac-4ddb-a6a9-55ee38a3e7ae
md" ### Quadratic approximation."

# ╔═╡ 7c8e8a61-b1c1-4f9a-addc-61a987c2c5d2
map_estimate = optimize(m4_1t, MAP())

# ╔═╡ a9d5118d-8fc7-4353-bb20-6b87ae0e9247
begin
	â, b̂, σ̂ = map_estimate.values
	(â = â, b̂ = b̂, σ̂ = σ̂)
end

# ╔═╡ f41d6486-7589-4b67-9eaf-92956622226a
let
	x = 30:0.1:65
	plot(x, â .+ b̂ .* x)
	scatter!(howell1.weight, howell1.height, leg=false)
end

# ╔═╡ 9f92d0a7-45ce-4148-b6a5-7ede9dc09234
md" ### Prediction."

# ╔═╡ 1cc80dac-1120-4520-b73f-0e31dc1d8721
begin
	x_test = [31, 40, 50, 60, 63]
	m_test = m4_1(x_test, fill(missing, length(x_test)))
	pred4_1t = predict(m_test, post4_1t)
	PRECIS(DataFrame(pred4_1t))
end

# ╔═╡ 159e07b6-dcd7-4ca7-9d3e-f5059b1b8803
md"##### Get the mean predicted values."

# ╔═╡ 85121bc0-447a-4b79-ac6c-ce91871ca7d8
begin
	pred_summary_df = DataFrame()
	pred_summary_df.parameters = names(pred4_1t)
	pred_summary_df.mean = mean(pred4_1t)[:, 2]
	pred_summary_df.std = vec(std(Array(group(pred4_1t, :heights)); dims = 1))
	pred_summary_df.lower = hpd(pred4_1t; alpha=0.11)[:, 2]
	pred_summary_df.upper = hpd(pred4_1t; alpha=0.11)[:, 3]
	pred_summary_df
end

# ╔═╡ 84b6ab83-7e65-477f-b87f-298d8b5eacea
let
	scatter(howell1.weight, howell1.height, lab="Observations",leg=:bottomright)
	scatter!(x_test, pred_summary_df.mean,
    	markersize=5, color=:red, lab="Predictions at x_test")
	x = range(minimum(howell1.weight), stop=maximum(howell1.weight), length=200)
	g(x) = mean(post4_1t[:a]) + mean(post4_1t[:b]) * x
	plot!(x, g.(x), lab="Regression line", leg=:bottomright)
end

# ╔═╡ 2d373995-f303-48b7-8cc3-1daf35d5ed36
md" ### StatisticalRethinking link() function."

# ╔═╡ 91a627e0-a46c-40f8-a5c2-2f3b25187074
@bind N Slider(1:1000, default=500)

# ╔═╡ a492bfeb-452c-4138-b5a2-6d2c4cdd114b
let
	xseq = 31:63
	μ = StatisticalRethinking.link(post4_1t_df, [:a, :b], xseq)
	μ = hcat(μ...)

	p = plot(; xlim=xseq, ylim=xseq, 
		xlab="weight", ylab="height", 
		title="Regression lines ($N) from draws",
		leg=false
	)
	for y ∈ first(eachrow(μ), N)
		plot!(p, xseq, y; c=:lightgrey, alpha=0.3)
	end
	p
end

# ╔═╡ ccf333ff-af94-439a-b41c-5be04c076007
md"""
!!! note

Click on "Live docs" and place cursor on link to see more help. 

Click little down arrow to the right to remove live docs again.
"""