using Distributions
using StatsPlots
using StatsBase
using KernelDensity

using StatisticalRethinking

# setting default attributes for plots
default(labels=false)

Code 3.1

Pr_Positive_Vampire = 0.95
Pr_Positive_Mortal = 0.01
Pr_Vampire = 0.001
tmp = Pr_Positive_Vampire * Pr_Vampire
Pr_Positive = tmp + Pr_Positive_Mortal * (1 - Pr_Vampire)
Pr_Vampire_Positive = tmp / Pr_Positive
Pr_Vampire_Positive

0.08683729433272395

Code 3.2

size = 1000
p_grid = range(0, 1; length=size)
prob_p = repeat([1.0], size);
prob_data = [pdf(Binomial(9, p), 6) for p in p_grid];
posterior = prob_data .* prob_p
posterior /= sum(posterior);

Code 3.3

samples_count = 10_000
cat = Categorical(posterior);
indices = rand(cat, samples_count)
samples = p_grid[indices];

Code 3.4

scatter(samples; alpha=0.2)

Code 3.5

density(samples)

Code 3.6

sum(posterior[p_grid .< 0.5])

0.17187458902022873

Code 3.7

sum(samples .< 0.5) / samples_count

0.1702

Code 3.8

sum(@. (samples > 0.5) & (samples < 0.75)) / samples_count

0.6034

Code 3.9

quantile(samples, 0.8)

0.7607607607607607

Code 3.10

quantile(samples, [0.1, 0.9])

2-element Vector{Float64}:
 0.44644644644644643
 0.8108108108108109

Code 3.11

size = 1000
p_grid = range(0, 1; length=size)
prob_p = repeat([1.0], size);
prob_data = [pdf(Binomial(3, p), 3) for p in p_grid];
posterior = prob_data .* prob_p
posterior /= sum(posterior)

samples_count = 10_000
cat = Categorical(posterior);
samples = p_grid[rand(cat, samples_count)];

Code 3.12

percentile(samples, [25, 75])

2-element Vector{Float64}:
 0.7087087087087087
 0.9329329329329329

Code 3.13

hpdi(samples, alpha=0.5)

2-element Vector{Float64}:
 0.8418418418418419
 0.998998998998999

Code 3.14

p_grid[argmax(posterior)]

1.0

Code 3.15

k = kde(samples, bandwidth=0.01)
k.x[argmax(k.density)]

0.9764975136347877

Code 3.16

mean(samples), median(samples)

(0.802118018018018, 0.8428428428428428)

Code 3.17

sum(@. posterior * abs(0.5 - p_grid))

0.3128751874998122

Code 3.18

loss = map(d -> sum(@. posterior * abs(d - p_grid)), p_grid);

Code 3.19

p_grid[argmin(loss)]

0.8408408408408409

Code 3.20

[pdf(Binomial(2, 0.7), n) for n ∈ 0:2]

3-element Vector{Float64}:
 0.09000000000000001
 0.41999999999999993
 0.4899999999999998

Code 3.21

rand(Binomial(2, 0.7))

1

Code 3.22

s = rand(Binomial(2, 0.7), 10)
println(s)

[2, 1, 2, 0, 1, 2, 1, 1, 1, 1]

Code 3.23

dummy_w = rand(Binomial(2, 0.7), 100_000);
proportions(dummy_w)  # or counts(dummy_w)/100000

3-element Vector{Float64}:
 0.08987
 0.41850000000000004
 0.49163000000000007

Code 3.24

dummy_w = rand(Binomial(9, 0.7), 100_000);
histogram(dummy_w; xlabel="dummy water count", ylabel="Frequency")

Code 3.25

w = rand(Binomial(9, 0.6), 10_000);

Code 3.26

w = [rand(Binomial(9, p)) for p in samples];