指数べき分布

・申請書の作成に午前中かかりきりになる．内容を変更するために，コピーしては貼り付けるだけなんだが，ワードの罫線はどうしてこう言うことを聞かないのか・・・．今回は関係ないけど，論文作成の時には，いつも”ページごと削除”というのができずに苦しんでいる．何とも情けない限りで，コンピュータに操られているとしか言い様がない．

・先日読んだRobledo-Arnuncio and Gil (2005)の論文の中で，「指数べき分布を使って花粉散布曲線を観察値にフィットさせる」というのがあるんだが，そもそも指数べき分布ってどんな感じだろうと思って，Rを使っての作図に取り組む．明らかにこんなことをやっている場合ではないのだが，忙しいときに限ってこんなことにはまってしまうわけで・・・．

・さて，指数べき分布は，a,bなどの2つのパラメータで決まる減少曲線である．Robledo-Arnuncioの論文にも式が書いてあるが，見慣れないΓ（ガンマと読む）という関数が出てくるのでエクセルでは手に負えない（・・・と思っていたら，Sさんは強引にΓ関数の値を算出する技を編み出している，すごい）．

・指数べき分布は，aとbの値によって，距離に対するスケールと形が異なるので，自由度が高く，「花粉散布や種子散布曲線の推定ではまず試みたまえ！」とClarkさんも述べている．平均値は，a×Γ（3/b)/Γ(2/b)で計算できるのだが，Rの場合，gamma()であっさりと計算してくれるので，式を書くのも楽だ，と分かる．

例えば，

a <- 20
b <- 0.7
exppower <- function(x) b/(2*pi*a^2*gamma(2/b))*exp(-(x/a)^b)
curve(exppower, 1, 1000, col="purple")

とでもすれば，a=20,b=0.7の時の指数べき分布曲線が紫色で描けるし，

avg.ab <- a*gamma(3/b)/(gamma(2/b))
print(a)
print(b)
print(avg.ab)

で平均値も求められる（もっとエクセレントな方法があるんだろうけど）．

aを動かしてみると，こんな感じである．b＝0.7に固定して，赤はa=10, 緑はa=20, 紫はa=30となっている．




今度はbを動かしてみる．a=20で固定．b=0.5が赤，b＝0.7が緑，b=0.9が紫．



・Robledo-Arnuncioの論文で載っていたa=24.08,b=0.67を代入して平均花粉散布距離を計算させると137.6である．たしか，こんな値だったわい，うまくできた・・・と思って論文を読み返すと，平均135.5mになっている．なぜであろうか・・・．小数点以下の問題かもしれないが妙に気になる．直接，本人にメールしてみるか，と思ったり・・・．

・ここまではいいのだが，最尤法で観察値にフィットさせるようにaとbを最適化させるのには，尤度方程式の理解が必要．未だに尤度という概念そのものが苦手である．まあとにかく，自分でやると勉強になるので，少しずつ進めることにする．