分布の裾あれこれ：Rでフライシュマンのべき乗変換・コーニッシュ・フィッシャー展開
キーワード
分布の歪み，高次積率，モーメント，尖度，歪度，テールリスク，ファットテール，VaR，期待ショートフォール，フライシュマンのべき乗変換，コーニッシュ・フィッシャー展開，Rスクリプト
背景：歪度と尖度ってなに
一般に外国語教育の研究者は，尖度や歪度のことをなめすぎだと思うの。縁あって，40年は続く学会誌の論文を全部読んでみたり（別の学会誌のは10年分は前に読んだ）したのだけど，尖度と歪度を報告している研究は，ほんとうに多くない。まあ，小標本の研究なんかで，頑張って標本の尖度や歪度を見ても…っていうのもわかるけど。でも…できれば報告しましょう。
そもそも，尖度と歪度は，高次積率（モーメント）とよばれる指標。
一次モーメントは，平均，二次モーメントは分散，三次モーメントが歪度，四次モーメントが尖度。三次以降を一般に高次という。
分散は，偏差の「二乗」の平均だから，二次。同じように，歪度は三乗，尖度は四乗する，っていう操作が入っている。
三次モーメントである歪度は，データを標準化したあとに，3乗し，その平均値をもとめた値。
四次モーメントである尖度は，データを標準化したあとに，4乗し，その平均値をもとめた値（または，そこから-3した数）
ちなみにn乗する時，nが奇数であれば左右対称な分布，偶数であれば，べき分布みたいになる。なぜかってのは，値の符号を考えればいい。
歪度は，大雑把だけども，平均では捉えられない位置関係（中心傾向）をあらわしていると思えば理解しやすい。平均よりも，データが平均よりも大きいところや小さいところに集まっていないか，みたいな。
尖度は標準偏差（分散）では捉えられないばらつき（散布度）をあらわしていると思えば理解しやすい。ピンとこない人もいるらしいのだけど，同じ標準偏差をもつデータでも，尖度の値が高ければ高いほどデータのばらつきが少なく，値が低ければ低いほど，データのばらつきが大きい。
もちろん，標準正規分布の歪度と尖度は0。
なんでそんなのが大事なの
金融とかの分野のひとには，お馴染みなのだろうけど，たとえばいきなり株が大暴落しないかとか，そういうのを予測して，管理して，意思決定するのが大事な場合もある。たとえば，損益の分布なんかは時系列にとると，儲けている時と損しているときなんか結構対称に出てきそう（なこともある）。そして，大損する確率も，大儲けする確率も少ない。大儲けするならいいけど，大損するリスクは考えもの。ちなみに，起きる確率は低いけども，それが起きると超困るようなものを，テールリスクなどというらしい。
これは，教育，たとえば処遇選択の意思決定にも似たようなことがいえる。
たとえば，まあほとんど起こらないだろうレベルの損失額の端っこを考えてもいい。正規分布するデータだったりしたら，予測区間の下限をもとめてもいい。これは金融とかリスク管理のひとは，VaR（value at risk）とよぶらしい。
一方，VaRは危ないとか，役に立たないっていう意見が主流になってきたんだそうだ。それで最近はVaRよりはむしろ，期待ショートフォール（expected shortfall）という指標を使うらしい。VaRとESについては，いつかちゃんと記事に書こう。
で，本題。
たとえば予測区間の下限は，平均と標準偏差からのみ計算されるのだけど，このとき，高次積率の情報を無視すると痛い目にあってしまう。たとえば，歪度が負だと左側の裾が重いので（ファットテールっぽい），標準偏差で予測される限界よりもはるかに，（まれに起こる）負の結果が起こりやすい。そして，尖度が低いと，標準偏差で予測される限界よりもはるかに（まれに起こる）負の結果が起こりやすい。
そりゃこまった。
外国語教育研究でもおなじ。
中心傾向だけで指導法の結果を議論すると，その裏で大部分のひとが泣いているかもしれない。
散布度（まで）だけで指導法の結果を議論すると，尖度や歪度の影響で正しい意思決定ができないかもしれない。
だから，尖度や歪度も大事にしましょうっていう話。記述統計量は等しく大事。
ごちゃごちゃいいからデータの例を出せや→フライシュマンのべき乗変換
はいはい。でも尖度や歪度を持ったデータを作るのって大変。いろいろなやり方があるけども，フライシュマンのべき乗変換っていうのをしたらオーソドックスでいいとおもわれ。フライシュマンのべき乗変換は，かならずってわけではないけど，ある程度任意の指定した尖度と歪度をもったデータを作ることができる。
細かいことは置いておいて，フライシュマン係数っていうのをまずはもとめなきゃならないのだけど，これは簡単じゃないので，Rのパッケージに頼る。"PoiNonNor"というパッケージで計算できる。
install.packages("PoisNonNor") #PoisNonNorパッケージ
library(PoisNonNor)
ここにある，Param.Fleishmanという関数にいれればフライシュマン係数を返してくれる。フライシュマン係数は4つセットででてくるね。
Param.fleishman(matrix(c(S,K),1,2) )
これ多変量対応なんだけど，とりあえず一変量だったら，Sに狙った歪度，Kに尖度を入れるといい。残念ながら収束しない組み合わせもある。
面倒なので，実際の変換などは関数にした。
rFPT<-function(n,M,SD,S,K){
SK<-matrix(c(S,K),1,2) #準備
Fc<-Param.fleishman(SK) #フライシュマン係数
x<-rnorm(n) #乱数の生成
x2<-Fc[1]+Fc[2]*x+Fc[3]*x^2+Fc[4]*x^3 #変換
x2<-x2*SD/sd(x2) #（不偏）標準偏差の調整
x2<-x2+(M-mean(x2)) #平均値の調整
x2 #結果
}
たとえば，1000個の，平均=0, 標準偏差 = 1，歪度 = 1，尖度 = 1に従う乱数なら，
rFPT(1000,0,1,1,1)->x
って感じ。ちょっと自分の都合だけど，平均と標準偏差は「かならず」指定した数値に一致するようにしているけども，歪度と尖度は一緒になるとは限らない。
できた変数のヒストグラムはこんな感じ。
うぃうぃ。
尖度や歪度の計算とか検定とか
尖度と歪度自体は簡単に計算できるけど，パッケージpsychの関数を使ってもいい。
install.packages("psych") #psychパッケージ
library(psych)skew(x)
kurtosi(x)
describe(x)
describe(x)$skew
describe(x)$kurtosis
せっかくなんで，専用のmomentsパッケージでもいい。
>install.packages("moments") #momentsskewness(x)
kurtosis(x) #これは3が基準！
非正規性の指標は，ジャック・ベラ統計量とか（JB）。
jarque.test(x)$statistic
検定は，なんでもいいけど，ダゴスティーノ検定とか，アンスコム・グリン検定とかをmomentsパッケージは推してた。
#ダゴスティーノ検定
agostino.test（x）#アンスコム・グリン検定
anscombe.test(x)
結局どうしたらいいの！
まあ，簡単にいえば，高次積率を加味した方法を使いましょうってこと。たとえば，コーニッシュ・フィッシャーの方法。
漸近展開のひとつ。漸近展開ってのは，超計算難しいやつを，厳密ではないけども，楽に計算する裏ワザやっちゃおうず，ってこと。
ここで使うコーニッシュ・フィッシャー展開は，尖度，歪度を与えて，累積正規分布の分位点を漸近的にもとめる。
まあ，簡単にいえば，分位点をもとめるとき，尖度と歪度を与えて，正規分布からのずれを補正しようっていうこと。
これ自体は計算が簡単なので関数作った。
CornFishPtile<-function(p,m,sd,s,k){
q<-qnorm(p) #累積標準正規分布から任意の分位点をもとめる
cf<-q+s*((q^2-1)/6)+k*((q^3-3*q)/24)-s^2*((2*q^3-5*q)/36) #コーニッシュ・フィッシャー展開
ptile<-m+cf*sd #平均値，標準偏差，推定値から値をもとめる
list("p"=p,"Standard_Normal"=q,"Cornish_Fisher_Estimate"=cf,"value"=ptile)
}
使い方は，平均 = 30, SD =5, 歪度 = -1，尖度 = 1の分布の95%点をもとめるとすると，
CornFishPtile(.95,30,5,-1,1)
だけでよい。
こんな感じで返す。
> CornFishPtile(.95,30,5,-1,1)
$percentile
[1] 0.95$Standard_Normal
[1] 1.644854$Cornish_Fisher_Estimate
[1] 1.321633$value
[1] 36.60816
$Standard_Normalが標準正規分布からもとめた値。$Cornish_Fisher_Estimateが補正した数値。結構違ってきちゃう。36.61くらいが95%点だろうってこと。
こういう仕組みを使ってVaRを計算することもあるらしい。modified VaRっていうらしいね。
でもこれやるくらいなら，はっきりいってブートストラップでもいいね…
分析例
たと えば，
rFPT(1000,3,.2,-.5,-.5)->x
こんなデータ。
>describe(x)
vars n mean sd median trimmed mad min max range skew kurtosis se
1 1 1000 3 0.2 3.02 3.01 0.23 2.42 3.3 0.87 -0.46 -0.56 0.01
>jarque.test(x)$statistic
JB
47.47136
これを，正規分布に則って，5%点を推定するとしたら，
mean(x)+qnorm(.05)*sd(x)
[1] 2.671029
これをコーニッシュ・フィッシャーを使うと，
> CornFishPtile(.05,3,.2,-.46,-.56)
$percentile
[1] 0.05$Standard_Normal
[1] -1.644854$Cornish_Fisher_Estimate
[1] -1.782939$value
[1] 2.643412
となるってこと。
まとめ
たとえば，歪んだデータのとき，標準偏差のみでリスクを評価するとリスクを過小評価してしまったりするかもしれない。中心傾向，そして散布度だけじゃなくて，歪みや尖りを考慮することも，強くもとめられると思う。特に，教育に関わるようなデータでは。
