尤度函数の形は単峰型で、山型の尤度函数の広がり方は$n$が大きくなればなるほど狭くなります。 サンプルをランダムに何回も生成し、尤度函数をプロットして視覚的に確認しましょう。
using Base64
displayfile(mime, file; tag="img") = open(file) do f
display("text/html", """<$tag src="data:$mime;base64,$(base64encode(f))">""")
end
using Plots
gr()
default(titlefontsize=10)
using ImplicitEquations
using Distributions
NaN
に設定)¶likelihood_func(p, q; x=[34, 32, 34]) =
if p + q ≤ 1
pdf(Multinomial(sum(x), [p, q, max(0, 1-p-q)]), x)
else
NaN
end
likelihood_func (generic function with 1 method)
function plot_random_linkelihood_func(n; p₀=0.3, q₀=0.3, α=0.05)
# サンプルを生成する確率分布
true_dist = Multinomial(n, [p₀, q₀, max(0, 1-p₀-q₀)])
# カイ二乗分布の相補累積分布関数
tau = cquantile(Chisq(2), α)
# サンプルXをランダムに生成
X = rand(true_dist)
k, l = X[1], X[2]
phat, qhat = k/n, l/n # 最尤法の解
MLE_str = "($(round(phat, digits=2)), $(round(qhat, digits=2)))"
p = q = range(0, 1, length=201)
# サンプルXから得られる尤度函数の値の計算
z = likelihood_func.(p', q; x=X)
P1 = plot(; title="n = $n, sample = ($k, $l, $(n-k-l))")
heatmap!(p, q, z; colorbar=false, aspectratio=1)
scatter!([p₀], [q₀]; color=:cyan, markerstrokewidth=0, markersize=3, label="true param")
scatter!([phat], [qhat]; color=:blue, markerstrokewidth=0, markersize=2.5, label=MLE_str)
plot!(xlim=(0,1), ylim=(0,1))
plot!(xlabel="p", ylabel="q")
P2 = surface(p, q, z; colorbar=false, gridalpha=0.5, grid_lw=0.5)
plot!(xlabel="p", ylabel="q")
plot(P1, P2; size=(800, 400))
end
plot_random_linkelihood_func (generic function with 1 method)
function gif_random_linkelihood_func(n; p₀=0.3, q₀=0.3, nframes=100, fps=2, fn="images/009_lik$n.gif")
anim = @animate for i in 1:nframes
plot_random_linkelihood_func(n; p₀=p₀, q₀=q₀, α=0.05)
end
gif(anim, fn, fps=fps)
displayfile("image/gif", fn)
end
gif_random_linkelihood_func (generic function with 1 method)
gif_random_linkelihood_func(10)
┌ Info: Saved animation to │ fn = C:\Users\4429s\codes\mathcodes\julia\statistics\images\009_lik10.gif └ @ Plots C:\Users\4429s\.julia\packages\Plots\6EMd6\src\animation.jl:104
gif_random_linkelihood_func(30)
┌ Info: Saved animation to │ fn = C:\Users\4429s\codes\mathcodes\julia\statistics\images\009_lik30.gif └ @ Plots C:\Users\4429s\.julia\packages\Plots\6EMd6\src\animation.jl:104
gif_random_linkelihood_func(100)
┌ Info: Saved animation to │ fn = C:\Users\4429s\codes\mathcodes\julia\statistics\images\009_lik100.gif └ @ Plots C:\Users\4429s\.julia\packages\Plots\6EMd6\src\animation.jl:104
gif_random_linkelihood_func(300)
┌ Info: Saved animation to │ fn = C:\Users\4429s\codes\mathcodes\julia\statistics\images\009_lik300.gif └ @ Plots C:\Users\4429s\.julia\packages\Plots\6EMd6\src\animation.jl:104
gif_random_linkelihood_func(1000)
┌ Info: Saved animation to │ fn = C:\Users\4429s\codes\mathcodes\julia\statistics\images\009_lik1000.gif └ @ Plots C:\Users\4429s\.julia\packages\Plots\6EMd6\src\animation.jl:104