草薙の研究ログ

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

Stanでカーブフィッティング:Menzerath-Altmann法則

もうしつこいのはわかっているけど,また,Menzerath-Altmann法則のはなし。この法則はよく知られた言語法則のひとつで,言語的構成要素の大きさは,その構成要素の大きさの関数だという観察で,具体的には,

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

という関数がよくフィットすることが知られている。

Menzerath-Altmann法則をRのnls関数で - 草薙の研究ログ


今日は,これをStanでやってみる。
MCMCによるカーブフィッティングは,このあたりの資料に詳しい。

Dugongs


ちょっとモデルを日本語で説明する。
観測データが,母数μとσをもつ正規分布に従う確率変数で,そのμは,

(a)μ = αx^{-β}e^{-λx}

という非線形関数に従っている。σは誤差。
α,β,λの事前分布はそれぞれ正規分布。σは,ガンマ分布に従うτから生成されているとする。
データ数が少ないから,それらしい事前分布を与えなければならない。具体的に,

  • α ~ Gaussian(0, 2)
  • β ~ Gaussian(0, 1)
  • λ ~ Gaussian(0, 1)
  • τ ~ Gamma(1/10000, 1/10000)

対象のデータはこんな感じ。

f:id:kusanagik:20180130201824p:plain

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サンプルの概観はこんな感じ。

f:id:kusanagik:20180130202055p:plain


これでMenzerath-Altmann法則はおしまい。長い前置きをして,やっとここまで来たという感じ。


以下のページを参考にさせていただいた。

Non-linear growth curves with Stan | mages' blog