Stanでカーブフィッティング:Menzerath-Altmann法則
もうしつこいのはわかっているけど,また,Menzerath-Altmann法則のはなし。この法則はよく知られた言語法則のひとつで,言語的構成要素の大きさは,その構成要素の大きさの関数だという観察で,具体的には,
(a)
という関数がよくフィットすることが知られている。
Menzerath-Altmann法則をRのnls関数で - 草薙の研究ログ
今日は,これをStanでやってみる。
MCMCによるカーブフィッティングは,このあたりの資料に詳しい。
ちょっとモデルを日本語で説明する。
観測データが,母数μとσをもつ正規分布に従う確率変数で,そのμは,
(a)
という非線形関数に従っている。σは誤差。
α,β,λの事前分布はそれぞれ正規分布。σは,ガンマ分布に従うτから生成されているとする。
データ数が少ないから,それらしい事前分布を与えなければならない。具体的に,
- α ~ Gaussian(0, 2)
- β ~ Gaussian(0, 1)
- λ ~ Gaussian(0, 1)
- τ ~ Gamma(1/10000, 1/10000)
対象のデータはこんな感じ。
MCMCは,反復回数 = 20000,焼却区間 = 5000,チェイン数 = 3くらいで。
これをやる。
#パッケージ library(rstan) library(HDInterval) #データ入力 dat<-list("N"=5, "x"=1:5, "Y"=c(3.1,2.53,2.29,2.12,2.09)) #プロット plot(s,m,xlab="Morpheme Length",ylab="Syllable Length",pch=20,cex=2) #Stan model <- " data { int<lower=0> N; real x[N]; real Y[N]; } parameters { real alpha; real beta; real lambda; real<lower=0> tau; } transformed parameters { real sigma; sigma = 1 / sqrt(tau); } model { real m[N]; for (i in 1:N) m[i] = alpha * pow(x[i],-beta) * pow(exp(1), -lambda * x[i]); Y ~ normal(m, sigma); alpha ~ normal(0, 2); beta ~ normal(0, 1); lambda ~ normal(0, 1); tau ~ gamma(.0001, .0001); } " #フィット fit <- stan(model_code = model, model_name = "MAL", iter=20000, warmup=5000, chains=3, data = dat) #結果を確認 fit #汚くてすみまそ par(mfrow=c(2,2)) plot(density(extract(fit)$alpha),type="n",main="alpha") polygon(density(extract(fit)$alpha),col="lightgreen") plot(density(extract(fit)$beta),type="n",main="beta") polygon(density(extract(fit)$beta),col="lightgreen") plot(density(extract(fit)$lambda),type="n",main="lambda") polygon(density(extract(fit)$lambda),col="lightgreen") plot(density(extract(fit)$sigma),type="n",main="sigma") polygon(density(extract(fit)$sigma),col="lightgreen") #最高密度区間で信用区間をチェック hdi(extract(fit))
それぞれの主要母数の事前分布に近似するMCMCサンプルの概観はこんな感じ。
これでMenzerath-Altmann法則はおしまい。長い前置きをして,やっとここまで来たという感じ。
以下のページを参考にさせていただいた。