草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

カーブフィッティング+ブートストラップ

また懲りずにMenzerath-Altmann法則の話。こんなデータがあって,

Morpheme Length(Syllables) Syllable Length(Phonemes)
1 3.10
2 2.53
3 2.29
4 2.12
5 2.09

Gabriel Altmann (1980). "Prolegomena to Menzerath's law". Glottometrika. 2: 1–10.


これがこの方程式にフィットするという観測的事実がある。

(a)y = ax^{-b}e^{-cx}

これを非線形最小二乗でやるのが普通なんだけど,今日はこれをブートストラップするという話。
MCMCやった後になんだけど,使うRパッケージは,nlstools。

#パッケージ読み込み
libarary(nlstools)

#データ入力
s<-1:5
m<-c(3.1,2.53,2.29,2.12,2.09)
dat<-data.frame(s,m)

#モデルをフィット
fit<-nls(m~a*s^(-b)*exp(1)^(-c*s),
start=list(a=3,b=0,c=0),
data=dat)

#ブートストラップ(B = 1,000)する
boot.fit<-nlsBoot(fit,1000)

#パーセンタイル法で誤差を見る
summary(boot.fit)

#適当にカーネルで可視化
plot(density(boot.fit$coefboot[,1]),main="a",type="n")
polygon(density(boot.fit$coefboot[,1]),col="lightblue")
plot(density(boot.fit$coefboot[,2]),main="b",type="n")
polygon(density(boot.fit$coefboot[,2]),col="pink")
plot(density(boot.fit$coefboot[,3]),main="c",type="n")
polygon(density(boot.fit$coefboot[,3]),col="lightgreen")

f:id:kusanagik:20180109205926p:plain

うっす。