読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

Tweedie分布のパラメータを推定する

@dichikaさんのブログ記事でTweedie分布の存在を知りました。Stanのメーリングリストでも「推定できないの?」という質問は過去にありましたが、多忙のBobさんからは「summing out(離散値をとるパラメータの和をとって消去)すればできるかもねー」という素っ気ない返答がありました。僕でもできたので備忘録として残します。

まずTweedie分布の紹介からします。

  • [1] Smyth, G. K. (1996) Regression modelling of quantity data with exact zeroes. Proceedings of the Second Australia-Japan Workshop on Stochastic Models in Engineering, Technology and Management. Technology Management Centre, University of Queensland, 572-580 (pdf file)
  • [2] Wikipedia (Web)
  • [3] SASのヘルプ (Web)
  • [4] Rの{twieedie}パッケージのdtweedie関数のヘルプ (Web)

密度関数の表し方やパラメータの書き方は様々あるようですが、この記事では[1]の記法を少しだけ変えた記法を使います。すなわち、 \lambda \alpha \beta([1]でいうところの \tauの逆数)を使って以下のように書きます( \lambda > 0,  \alpha > 0,  \beta > 0を満たします)。

1.png

 Y = 0の場合は確率  e^{- \lambda } で生成されます。 Y > 0の場合は複合ポアソン-ガンマ分布(compound Poisson–gamma distribution)から生成されます。ガンマ分布には再生性がありますので、以下のように書くこともできます(参考記事)。

2.png

なぜ \tauではなく \betaに変えているかというと、ガンマ分布はscaleパラメータで指定する方法とrateパラメータで指定する方法があって、StanやBUGSではrateパラメータを使って指定し、論文ではscaleパラメータを使って指定しているからです(Wikipedia参照)。ちなみにRのdgamma関数はどちらでも指定可能です。

Rのtweedieパッケージを使う場合は、Tweedie分布のパラメータは \mu \phi \xi([1]でいうところの \theta)で指定します。Tweedie分布の平均の理論値は \mu、分散の理論値は \phi \mu^{\theta}になります。 \mu > 0,  \phi > 0,  1 \lt \theta \lt 2を満たします。なお、 \mu \phi \theta \lambda \alpha \betaの間には以下のような関係があります。

3.png

4.png

パラメータを色々試して密度関数も見ておきましょう。

tweedie.png

上の図では見やすさのため Y = 0のところだけ棒グラフにしています。Zero-inflatedな分布の一つです。特に \thetaが1に近い場合は、ポアソン分布の離散性と分散の小さいガンマ分布の複合が効いて変態的な分布となっています。ちなみにパラメータの値は以下の通りです。

muphithetalambdaalphabeta
0.30.21.011.5399506.06
10.21.015.0599500
30.21.0114.9999494.54
0.30.51.010.6199202.42
10.51.012.0299200
30.51.015.9999197.81
0.311.010.3199101.21
111.011.0199100
311.0139998.91
0.331.010.19933.74
131.010.349933.33
331.0119932.97
0.30.21.051.6819106.2
10.21.055.2619100
30.21.0514.951994.66
0.30.51.050.671942.48
10.51.052.111940
30.51.055.981937.86
0.311.050.341921.24
111.051.051920
311.052.991918.93
0.331.050.11197.08
131.050.35196.67
331.051196.31
0.30.21.33.082.3323.92
10.21.37.142.3316.67
30.21.315.412.3311.99
0.30.51.31.232.339.57
10.51.32.862.336.67
30.51.36.162.334.79
0.311.30.622.334.78
111.31.432.333.33
311.33.082.332.4
0.331.30.212.331.59
131.30.482.331.11
331.31.032.330.8
0.30.21.944.330.1116.42
10.21.9500.115.56
30.21.955.810.112.07
0.30.51.917.730.116.57
10.51.9200.112.22
30.51.922.320.110.83
0.311.98.870.113.28
111.9100.111.11
311.911.160.110.41
0.331.92.960.111.09
131.93.330.110.37
331.93.720.110.14

Stanで実装する際は、複合ポアソン-ガンマ分布の Mが離散値をとるパラメータなので、和をとって変数消去する必要があります。変数消去はたしかに面倒ですが、(Stanに限らず)推定が高速になるメリットがあります。

今回あてはめたいデータの分布と上記の密度関数の図を見比べると、Tweedie分布のパラメータthetaの値は1.3以下になりそうだと判断できたとします。すると、パラメータlambdaは15程度、すなわちMとしては1~30までを見ておけば十分によく近似できそうです。

以上よりStanで実装した例は以下になります。

  • 24~25行目:  Y = 0の場合は確率  e^{- \lambda } で生成されます。Stanではtarget記法を使って対数確率をtargetに加えています。ベルヌーイ分布を使うこともできるでしょう。
  • 27~30行目:  Y > 0の場合はsumming outしてます。説明をブログでするのは面倒なので、log_sum_exp関数については、Stan 2.9.0のマニュアルの「10. Finite Mixtures」章の「10.3. Log Sum of Exponentials: Linear Sums on the Log Scale」節を参考にしてもらえればと思います。

このモデルではmu,phi,thetaを推定してlabmda,alpha,betaを算出しています。逆に、labmda,alpha,betaparametersブロックに置いて推定し、mu,phi,thetagenerated quantitiesブロックに置いて算出するモデルも試しましたが、alpha,betaがとても大きな値になる可能性があるために、これらに対する弱情報事前分布が設定できず、推定がとても遅くなりました。

データを作成して、推定を実行するRコードは以下になります。ここでは2種類のパラメータの組み合わせから乱数を生成しました。

推定結果は以下になります。

data1

meanse_meansdX2.5.X25.X50.X75.X97.5.n_effRhat
mu0.970.000.070.840.920.971.021.1118971.00
phi1.000.000.100.830.941.001.071.2111121.00
theta1.310.000.041.251.291.311.341.3915701.00
lambda1.430.000.131.201.351.431.521.7020361.00
alpha2.230.010.381.571.962.202.463.0615571.00
beta3.300.020.652.252.843.243.674.8313421.00
lp__-285.970.041.28-289.35-286.54-285.61-285.05-284.5612281.00

data2

meanse_meansdX2.5.X25.X50.X75.X97.5.n_effRhat
mu3.020.000.062.912.993.023.063.1315801.00
phi1.000.000.001.001.001.001.001.0040001.00
theta1.010.000.001.011.011.011.011.0116991.00
lambda3.020.000.062.912.983.023.063.1316131.00
alpha94.440.114.6185.5691.4594.3397.51103.6917181.00
beta94.400.114.6285.5191.4094.2997.50103.5617211.00
lp__-1581.310.031.23-1584.42-1581.86-1580.98-1580.43-1579.9413231.00

回帰にも少しいじれば使えると思います。