遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

安定分布は４つのパラメタで描く

R 安定分布 分布
  • 2012/12/06にうまくfBasicsパッケージの安定分布関数が動かないので…
  • stabledistパッケージのdstable()関数などを使う方がよいかも…
  • 安定分布は"skewness and heavy tails"な分布を包含する分布群
  • 安定分布に関する基本的な論文を書き、ホームページを開設しているJohn P. Nolan氏のホームページ
  • 経済分野での活用が活発な模様
  • 特性関数で一意に定義できる
  • Rには安定分布のためのいわゆる分布関数("dstable()","pstable()","qstable()","rstable()"がある。fBasicsパッケージに入っている
  • 安定分布は(安定分布の特性関数は)４つのパラメタで定義される
    • \alpha,0 < \alpha \le 2：分布の厚みを定める。小さいと裾が広い
    • \beta,-1 \le \beta \le 1：歪度指数、あるいは非対称パラメータ。0のときに左右対称
    • \gamma, \gamma >0：規模のパラメタ。拡大縮小する。
    • \delta：平行移動のパラメタ。
install.packages("stabledist")
library(stabledist)
x<-seq(from=-10,to=10,by=0.1)
plot(x,dstable(x,alpha=1,beta=0.99,gamma=1,delta=2),type="l")
  • 安定分布の特性関数は次のように表される
    • exp(i\delta z - \gamma |z|^{\alpha} (1+i \beta sgn(z) w(z,\alpha)))
    • w(z, \alpha)=\tan (\frac{\pi \alpha}{2}), \alpha \not=1
    • w(z, \alpha)=\frac{\pi}{2} \log(|z|), \alpha = 1
  • この表現から、ある意味で余計な部分を取り去る
    • 平行移動・拡縮は不要、左右非対称も取り去ると
    • exp(\gamma |z|^{\alpha}
    • レヴィ分布は\alpha=0.5で、\beta=1と左右非対称性を最大化したもの
  • ここで、『安定分布は、ある分布に従う、複数の独立な確率変数の１次線形結合した(足し合わせた)ときに、部分を構成している分布になるという特徴を持って定義されている』
  • また、独立な確率変数の１次線形結合にあっては、特性関数の掛け算で結合後の特性関数が得られるから、\alpha, \betaが等しい特性関数２個の足し合わせの特性関数は以下のようになって、『特性関数の掛け合わせが、\deltaと\gammaの足し合わせの部分だけが変化して、\alphaと\betaは変わらない。つまり、\alpha,\betaが確率変数の大事な特徴を決めているときには、「足した後もその大事な特徴が維持される」ことを意味している。\gamma,\deltaは拡縮・平行移動のパラメタであって、「形の特徴」には影響せず、\alpha,\betaが「形の特徴」と決めていると、上で書いてあるのを思い出せば、腑に落ちる
    • exp(i(\delta_1+\delta_2) z - (\gamma_1 + \gamma_2) |z|^2 (1+i \beta sgn(z) w(z,\alpha)))
  • 安定分布なら、足しても分布の形が変わらないけれど、安定分布でない分布(特性関数が安定分布の特性関数の形をしていない関数(たとえば指数関数)では、足すと分布の形が変わる様子を試すソースが以下
a<-runif(1)*2
b<-(runif(1)-0.5)*2
#a<-1
#b<-0.99
N<-100
X1<-rstable(N,alpha=a,beta=b)
X2<-rstable(N,alpha=a,beta=b)
X12<-X1+X2

Y1<-rexp(N,0.1)
Y2<-rexp(N,0.1)
Y12<-Y1+Y2
layout(matrix(c(1,2,3,4,5,6),2,3,byrow=TRUE))
plot(sort(X1))
plot(sort(X2))
plot(sort(X12))

plot(sort(Y1))
plot(sort(Y2))
plot(sort(Y12))

# 偏った分布として安定分布をもとの分布としてとる
#x <- rexp(100000)
install.packages("stabledist")
library(stabledist)

x <- rstable(100000,alpha=0.6,beta=0.9,gamma=1,delta=0)
# その分布を見てみる
hist(x)
# サンプリングする数を1,2,4,...と変化させ、
# 1000回ずつやって、そのサンプル平均の分布を見てみよう
n.iter <- 1000
n.sample <- 2^(0:15)
# 1000回ずつの平均を格納する
mean.stocker <- matrix(0,length(n.sample),n.iter)
#par(ask = TRUE)
for(i in seq(n.sample)){
# n.sample[i]個ずつn.iter回サンプルするのを一度にやって
# それを行列にする
	tmp.X <- matrix(sample(x,n.iter*n.sample[i],replace=TRUE),nrow = 
n.iter)
# 行ごとに(試行ごとに)平均をとる
	means <- apply(tmp.X,1,mean)
	hist(means)
	mean.stocker[i,] <-sort(means)
}
# ソートしてプロットする
# サンプル数が増えるほど、ほとんど同じ平均値が得られることが
# グラフが横に寝ることでわかる
matplot(t(mean.stocker),type="l")
# ゆっくりと、平均値のヒストグラムを描くと
# たしかに、だんだん正規分布化する
par(ask =TRUE)
for(i in seq(n.sample)){
	title <- c("N sampled is ",n.sample[i])
	hist(mean.stocker[i,],main=title)
}
par(ask = FALSE)