ニュートン法（ガンマ分布の最尤推定）

MLE_Gamma.jl

using Distributions
using SpecialFunctions
using Random
using Plots, StatsPlots

function Mstep(d::Gamma, y)
    shp, scl = params(d)
    rho = log(shp)
    meanlogy = mean(log,y)
    #f(x,shp) = -shp*log(scl) + (shp-1)*log(x) - (x/scl) -loggamma(shp)
    g = shp*log(scl) - shp*meanlogy + digamma(shp)*shp
    H = shp*log(scl) - shp*meanlogy + trigamma(shp)*shp^2 + digamma(shp)*shp
    rho -= g/H
    b = mean(y)/exp(rho)
    return Gamma(exp(rho), b)
end

function mynewton(d, y, iter)
    p = zeros(iter, 2)
    for i in 1:iter
    d = Mstep(d, y)
    p[i,:] .= log.(params(d))
    end
    return p
end

y = rand(MersenneTwister(12345), Gamma(2.0,2.0), 100)
theta = mynewton(Gamma(1.0, 1.0), y, 100)
maximum(theta[:,1])
minimum(theta[:,1])
maximum(theta[:,2])
minimum(theta[:,2])
av = -1.5:0.05:3.5
bv =  -1.5:0.05:3.5
Xv = repeat(reshape(av, 1, :), length(bv), 1)
Yv = repeat(bv, 1, length(av))
Zv = map((a,b) -> exp(mean(y -> logpdf(Gamma(exp(a),exp(b)),y), y)), Xv, Yv)

anim = @animate for i in 1:100
    contour(av, bv, Zv, linewidth=1, cb=false)
    scatter!([theta[i,1]],[theta[i,2]],ms=6,color="black",legend=false)
end

gif(anim, "./Desktop/anim_fps15.gif", fps = 15)