草薙の研究ログ

英語の先生をやってます。

当てはめた分布と推定した母数から平均,分散,歪度,尖度をもとめる

なんかある論文で,分布が強く歪んでいることが理論的に明白だった変数に,ガンマ分布か対数正規分布かなにかを当てはめて,その推定母数と適合度指標のみを(きつい紙幅の関係もあって)報告したときに,「標本の記述統計を報告しないとはけしからん」というようなことを査読者さまにご教示されたことがある。ま,ごもっともであるけど,紙幅の関係もあって,ううむ,標本自体がそんな重要なのかな…みたいに思いつつも,どうしたんだったかな。結構前の話だ。

でも,もちろん標本の記述統計量はわからんけど,この当てはめた分布の母数の推定値から,理論的には,そしてその分布モデルがしっかりと当てはまっているなら,標本のおおよその記述統計であれば,あとからでもとめることもできる。たとえばガンマ分布はかなり簡単で,こんな自作関数でも作ればいい。

gammadescriptive<-function(shape,scale){
	m<-shape*scale
	v<-shape*(scale^2)
	s<-2/sqrt(shape)
	k<-6/shape
	result<-list("mean"=m,"variance"=v,"skew"=s,"kurtosis"=k)
	result
}

たとえば,この関数に,形状母数3,尺度母数(scale)4を入れてやる。

gammadescriptive(3,4)

$mean
[1] 12

$variance
[1] 48

$skew
[1] 1.154701

$kurtosis
[1] 2

するとこんな感じになる。もちろん,「あくまでも理論値」なんだけど。

こんな感じで実験してみよう。

#パッケージの準備
library(psych)
library(fitdistrplus)

#例の作成
set.seed(1)
dat<-rgamma(1000,shape=3,scale=4)

#理論値
gammadescriptive(3,4)

$mean
[1] 12

$variance
[1] 48

$skew
[1] 1.154701

$kurtosis
[1] 2


#この例の記述統計を見てみる
describe(dat)
   vars    n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
X1    1 1000 11.95 7.09   10.4   11.19 6.59 0.74 50.2 49.46 1.12     1.86 0.22

#大体近い
#また,この乱数にガンマ分布をフィットさせる
model<-fitdist(dat,"gamma")
model
Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters:
       estimate Std. Error
shape 2.8585896 0.12111673
rate  0.2392977 0.01108254

#ちっ,これrateで出してきやがる,でもscale=1/rate
#でも,大体形状母数3,尺度母数4くらい

sh<-model$estimate[[1]]
sc<-1/model$estimate[[2]]

#またこの母数の理論値を見てみる
gammadescriptive(sh,sc)
$mean
[1] 11.94575

$variance
[1] 49.92004

$skew
[1] 1.182917

$kurtosis
[1] 2.098937


もちろんだけど,おおよそ標本の記述統計にも合っている。
だから,ある変数に狙った分布がよく当てはまっているといえるのであれば,平均,分散,歪度,尖度などは報告しなくても,「後からおおよその値は計算できる」こともあることは事実。ただ,この記述統計4つを報告したところで,この経験分布はガンマ分布に近いとか適合度はどれくらいだとか,そういうのはわからないか,または技術的により困難だってこと。

もちろん,標本の記述統計は要りません,って話じゃないけど。