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

## モチベーション
1次データが手に入らずおおまかな度数分布しかないとき、限られたデータ（度数分布／分位数）からサンプルのGamma分布をTuringで調べたい。
![gamma-question.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/ba4d7fc7-513c-6a58-2c51-4a692c6bc39f.png)

## Gamma分布！

https://qiita.com/kimseok1973/items/1bff01083e0f7e2ccde7

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

```julia
    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の逆関数。

```Julia
    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データを次のように準備する。
```Julia
    q4 = quantile(Gamma(3,1.5),0.4) #q4= 3.427615356005071
    q2 = quantile(Gamma(3,1.5),0.2) #q2= 2.302566303966965
```
モデルを次のように定義して実行する。
```Julia
    @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           | 

![marathon-3graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/f68009ab-f30c-10b5-c90c-1a0885984369.png)