Embedded Files
# Chapter Appendix A
#import Pkg
#Pkg.add("Plots")
#Pkg.add("Optim")
#Pkg.add("StatsKit")
#Pkg.add("GLM")
#Pkg.add("DataFrames")
#Pkg.add("StatsPlots")
#Pkg.add("StatFiles")
#Pkg.add("NLSolversBase")
#Pkg.add("Distributions")
using Random
using Statistics
using Plots
using Optim
using NLSolversBase
using StatsKit
using CSV
using GLM
using DataFrames
using LinearAlgebra
using StatsPlots
using StatFiles
using Distributions
# A.1
# A.2
# A.2.1
# A.2.2
rng = MersenneTwister(123456789)
N = 10
μ = -3
σ = 5
x = randn!(rng, zeros(N)).*σ .+ μ
dens_x = kde(x)
mean(x)
std(x)
plot(dens_x)
plot!(kde(randn!(rng, zeros(10000)).*σ .+ μ))
# A.2.3
rng = MersenneTwister(123456789)
M = 500
sample_est = zeros(M, N)
for m = 1:M
sample_est[m,1:end] = randn!(rng, zeros(N)).*σ .+ μ
end
means = zeros(M)
for m = 1:M
means[m] = mean(sample_est[m,1:end])
end
histogram(means)
vline!(fill(-3,M))
# A.2.4
rng = MersenneTwister(123456789)
N = 100000
x = randn!(rng, zeros(N)).*σ .+ μ
mean(x)
# A.2.5
rng = MersenneTwister(123456789)
M = 1000
sample_means = zeros(M,4)
for m = 1:M
sample_means[m,1] = mean(randn!(rng, zeros(10)).*σ .+ μ)
sample_means[m,2] = mean(randn!(rng, zeros(100)).*σ .+ μ)
sample_means[m,3] = mean(randn!(rng, zeros(1000)).*σ .+ μ)
sample_means[m,4] = mean(randn!(rng, zeros(10000)).*σ .+ μ)
end
plot(kde((sample_means[1:end,1] .- μ).*(√9/σ)))
plot!(kde((sample_means[1:end,2] .- μ).*(√99/σ)))
plot!(kde((sample_means[1:end,3] .- μ).*(√999/σ)))
plot!(kde((sample_means[1:end,4] .- μ).*(√9999/σ)))
# Note here I have adjusted for degrees of freedom.
# A.2.6
rng = MersenneTwister(123456789)
M = 5
N = 10
sample_means = zeros(M,2)
for m = 1:M
x₁ = randn!(rng, zeros(N)).*σ .+ μ
sample_means[m,1:end] = [mean(x₁); std(x₁)/(√(N-1))]
end
plot(kde((randn!(rng, zeros(1000)) .+ sample_means[1,1]).*sample_means[1,2]))
plot!(kde((randn!(rng, zeros(1000)) .+ sample_means[2,1]).*sample_means[2,2]))
plot!(kde((randn!(rng, zeros(1000)) .+ sample_means[3,1]).*sample_means[3,2]))
plot!(kde((randn!(rng, zeros(1000)) .+ sample_means[4,1]).*sample_means[4,2]))
plot!(kde((randn!(rng, zeros(1000)) .+ sample_means[5,1]).*sample_means[5,2]))
# A.2.8
rng = MersenneTwister(123456789)
Random.seed!(123456789)
N = 10
x₁ = randn!(rng, zeros(N)).*σ .+ μ
M = 1000
bs_mean = zeros(M)
for m = 1:M
bs_mean[m] = mean(rand(x₁, N))
end
plot(kde(bs_mean))
plot!(kde(randn!(rng, zeros(100000)).*(σ/√(N-1)) .+ μ))
# A.2.9
rng = MersenneTwister(123456789)
α = 0.05
cr = std(x).*quantile(randn!(rng, zeros(10000)), α)./√9
# Again adjusting for degrees of freedom.
plot(kde(randn!(rng, zeros(10000)).*(std(x)/√9)))
vline!(fill(cr,10))
vline!(fill(mean(x),10))
mean(randn!(rng, zeros(10000)) .< √9*mean(x)/std(x))
# A.3
# A.3.1
# A.3.2
# A.3.3
logϕ = function(z)
-z.^2 - √π
end
log_norm = function(x, μ, σ)
z = (x .- μ)./σ
sum(logϕ.(z) .- log(σ))
end
Θ = ((1:100).*15)./100 .- 10
g = (1:100)./100
h = []
for i = 1:100
h = [h; log_norm(x₁, Θ[i], 5)]
end
γ = exp.(h .+ log.(g))./sum(exp.(h .+ log.(g)))
plot(Θ, γ)
vline!(fill(-3,100))
vline!(fill(mean(x₁), 100))
# A.4
# A.4.1
# A.4.2
# A.4.3
Θ = (maximum(sample_est) - minimum(sample_est)).*((1:1000)./1000) .+ minimum(sample_est)
H = zeros(500,1000)
for i = 1:500
x = sample_est[i,1:end]
for j = 1:1000
H[i,j] = log_norm(x, Θ[j], 5)
end
end
g₀ = fill(1/1000, 1000)
log_G0 = log.(zeros(500, 1000) .+ 1/1000)
γ₀ = exp.(H + log_G0)
for i = 1:500
γ₀[i,1:end] = γ₀[i,1:end]./sum(γ₀[i,1:end])
end
g₁ = []
for i = 1:1000
g₁ = [g₁; sum(γ₀[1:end,i])]
end
g₁ = g₁./sum(g₁)
log_G1 = log.(g₁)
for i = 2:500
log_G1 = hcat(log_G1, log.(g₁))
end
γ₁ = exp.(H + log_G1')
for i = 1:500
γ₁[i,1:end] = γ₁[i,1:end]./sum(γ₁[i,1:end])
end
g₂ = []
for i = 1:1000
g₂ = [g₂; sum(γ₁[1:end,i])]
end
g₂ = g₂./sum(g₂)
plot(Θ, g₀)
plot!(Θ, g₁)
plot!(Θ, g₂)
ϵ = 1e-3 # small number
g_old = g₁
g_new = g₂
diff = sum(abs.(g_old - g_new))
while diff > ϵ
g_old = g_new
log_G = log.(g_old)
for i = 2:500
log_G = hcat(log_G, log.(g_old))
end
γ = exp.(H + log_G')
for i = 1:500
γ[i,1:end] = γ[i,1:end]./sum(γ[i,1:end])
end
g_new = []
for i = 1:1000
g_new = [g_new; sum(γ[1:end,i])]
end
g_new = g_new./sum(g_new)
diff = sum(abs.(g_old - g_new))
println(diff)
end
μ_est = []
for i = 1:500
μ_est = [μ_est; mean(sample_est[i])]
end
J = 100000
g_dist = []
for i = 1:1000
if round(J*g_new[i]) > 0
g_dist = [g_dist; fill(Θ[i], Int(round(J*g_new[i])))]
end
end
plot(density(g_dist))
vline!(fill(-3, 100))
plot!(density!(μ_est))
# A.4.4
h_x1 = []
for i = 1:1000
h_x1 = [h_x1; log_norm(x₁, Θ[i], 5)]
end
γ_x1 = (exp.(h_x1) + g_new)./sum(exp.(h_x1) + g_new)
J = 100000
g_dist_x1 = []
for i = 1:1000
g_dist_x1 = [g_dist_x1; fill(Θ[i], Int(round(J*γ_x1[i])))]
end
plot(density(g_dist_x1))
vline!(fill(-3, 100))
vline!(fill(mean(x₁), 100))
# A.5
# A.5.1
# A.5.2
log_bin = function(N, p̂, p)
K = Int(round(p̂*N))
return log(binomial(big(N), big(K))) + K*log(p) + (N-K)*log(1 - p)
end
Θ = (1:999)./1000
log_g = fill(log(1/1000), 999)
h_jp = []
for i = 1:999
h_jp = [h_jp; log_bin(3, 1, Θ[i])]
end
γ_jp = exp.(h_jp + log_g)./sum(exp.(h_jp + log_g))
J = 100000
γ_dist_jp = []
for i = 1:999
γ_dist_jp = [γ_dist_jp; fill(i, Int(round(J*γ_jp[i])))]
end
# mean
γ_jp'*Θ
plot(density(γ_dist_jp))
vline!(fill(342,100))
xlabel!("batting average")
sum(γ_jp[342:999])
# A.5.3
M = 1000
path = "ENTER PATH HERE"
x = CSV.read(string(path,"/bat_ave.csv"))
# data set created from Lahman's data
H = zeros(M, 999)
for m = 1:M
for i = 1:999
H[m,i] = log_bin(x.AB[m], x.AVE[m], Θ[i])
end
#println(m)
end
log_G = fill(log(1/1000), M, 999)
γ₀ = exp.(H + log_G)
for m = 1:M
γ₀[m,1:end] = γ₀[m,1:end]./sum(γ₀[m,1:end])
end
g₁ = []
for i = 1:999
g₁ = [g₁; sum(γ₀[1:end,i])]
end
g₁ = g₁/sum(g₁)
ϵ = 0.001
g_old = exp.(log_g)
g_new = g₁
diff = sum(abs.(g_new - g_old))
while diff > ϵ
g_old = g_new
log_G = log.(g_old)
for i = 2:M
log_G = hcat(log_G, log.(g_old))
end
Γ = exp.(H + log_G')
for m = 1:M
Γ[m,1:end] = Γ[m,1:end]./sum(Γ[m, 1:end])
end
g_new = []
for i = 1:999
g_new = [g_new; sum(Γ[1:end,i])]
end
g_new = g_new./sum(g_new)
diff = sum(abs.(g_old - g_new))
println(diff)
end
# A.5.4
γ_jp2 = exp.(h_jp + log.(g_new))
γ_jp2 = γ_jp2./sum(γ_jp2)
J = 100000
g_dist_mlb = []
for i = 1:999
g_dist_mlb = [g_dist_mlb; fill(Θ[i], Int(round(J*g_new[i])))]
end
γ_dist_jp2 = []
for i = 1:999
γ_dist_jp2 = [γ_dist_jp2; fill(Θ[i], Int(round(J*γ_jp2[i])))]
end
# mean
γ_jp2'*Θ
sum(γ_jp2[342:999])
plot(density(g_dist_mlb), label="Prior")
plot!(density!(γ_dist_jp2), label="JPs Posterior")
vline!(fill(0.342, 100))
xlabel!("batting average")
# A.6
# A.6.1
# A.6.2
# A.6.3
# A.6.4
# A.7
Google Sites
Report abuse