Turingで度数分布からGamma分布のパラメータを調べる

モチベーション

1次データが手に入らずおおまかな度数分布しかないとき、限られたデータ（度数分布／分位数）からサンプルのGamma分布をTuringで調べたい。

Gamma分布！

この記事の続きです。マラソンランナーの到着時間を標準分布で推定しましたが、それぞれ独立した事象（ランナー）の到着時間を表すガンマー分布のほうが適したモデルかもしれない。そこで度数データで推定したいのですが、実はJuliaのgammacdf関数が自動微分をサポートしないので、前項で示した分位数から推定する方法ではエラーになってしまいます。

    MethodError: no method matching __gamma_inc_inv(::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 4}, ::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 4}, ::Bool)
    Closest candidates are:
      __gamma_inc_inv(::T, ::T, ::Bool) where T<:Union{Float16, Float32} at C:\Users\xxxxxx\.julia\packages\SpecialFunctions\hefUc\src\gamma_inc.jl:1089
      __gamma_inc_inv(::Float64, ::Float64, ::Bool) at C:\Users\xxxxxx\.julia\packages\SpecialFunctions\hefUc\src\gamma_inc.jl:1012
    
    Stacktrace:....省略

どうやらSpecialFunctionsパッケージの関数__gamma_inc_invが自動微分（ForwardDiff）のstructを認識できないのでこれを強引に食わせることで解決をはかります。gamma_inc_invは不完全ガンマ函数gamma_incの逆関数。

    using ForwardDiff,StatsFuns,SpecialFunctions
    function SpecialFunctions.__gamma_inc_inv(a::ForwardDiff.Dual{T}, 
            minpq::ForwardDiff.Dual{T}, pcase::Bool) where {T}
        v_a= ForwardDiff.value(a)
        v_pq = ForwardDiff.value(minpq)
        SpecialFunctions.__gamma_inc_inv(v_a,v_pq,pcase)
    end

Toyデータを次のように準備する。

    q4 = quantile(Gamma(3,1.5),0.4) #q4= 3.427615356005071
    q2 = quantile(Gamma(3,1.5),0.2) #q2= 2.302566303966965

モデルを次のように定義して実行する。

    @model function gdemo()
        m ~ InverseGamma(2,3)
        s ~ InverseGamma(2,3)
        dist = Gamma(m,s)
        
        y4 ~ Normal(quantile(dist,0.4),0.1)
        y2 ~ Normal(quantile(dist,0.2),0.1)
        
        y4 ~ Normal(q4,0.1) #q4 is a constance , 3.4276
        y2 ~ Normal(q2,0.1) #q2 is a constance , 2.3025
    end   
    
    c = sample(gdemo(),NUTS(),6000)

結果は次の通りで正解に近似した結果が得られた。

Summary Statistics
      parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
          Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯
    
               m    3.0844    0.8489     0.0110    0.0841   46.3321    1.0234      ⋯
               s    1.6008    0.5086     0.0066    0.0485   52.3757    1.0242      ⋯
                                                                    1 column omitted
    Quantiles
      parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
          Symbol   Float64   Float64   Float64   Float64   Float64 
    
               m    1.9430    2.4611    2.8759    3.5576    5.0975
               s    0.7607    1.2086    1.5787    1.9368    2.6341

quantileはうまくいくがcdf関数はうまくいかない。

gammacdf関数が呼び出すgamma_incのAuto Difference対応は難しそうでお手上げ。しかし知っている情報が階級数と分位数だけなら、cdfでもquantileでも構わないのでとりあえず良しとする。

マラソンのランナー分布

前回の課題であった北海道マラソンのランナー到着時間分布をGamma分布で推定した結果を掲載する。ガンマ分布を採用した場合95％の確信区間として211分～317分と推定できる。

	累積度数分布	標準分布	ガンマ分布
1	0.001820	0.002075	0.001080
2	0.003389	0.004362	0.003015
3	0.008913	0.008685	0.007355
4	0.017324	0.016384	0.015944
5	0.035652	0.029310	0.031169
6	0.054105	0.049763	0.055634
7	0.085865	0.080268	0.091634
8	0.130304	0.123147	0.140568
9	0.179450	0.179960	0.202470
10	0.238137	0.250907	0.275820
11	0.319232	0.334413	0.357692
12	0.844652	0.849784	0.815851
13	1.0	0.999963	0.998646