草薙の研究ログ

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

brmsパッケージでカーブフィッティング

計量言語学のMenzerath-Altmann法則をbrmsで。

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

s<-1:5
m<-c(3.1,2.53,2.29,2.12,2.09)
dat <- data.frame(s, m)
my.prior <- prior(normal(0,3), nlpar = a) + prior(normal(0,1), nlpar = b) + prior(normal(0,1), nlpar = c)
fit <- brm(bf(m ~ (a*s^-b)*exp(1)^(-c*s), a + b + c ~ 1, nl = TRUE),
	data = dat,
	prior = my.prior,
	save_all_pars=T)
summary(fit)


f:id:kusanagik:20180914171856p:plain


f:id:kusanagik:20180914172534p:plain

これ,データポイントが少ないから,事前分布の設定がすげえ問題だよね。
適当に,

  • a ~ Normal(0,3)
  • b ~ Normal(0,1)
  • c ~ Normal(0,1)

とかにした。実際に-1とか1とかの値を取ることがほとんどないbcはともかく,このaの設定は恣意的といわざるをえないけど。
ちなみに,より簡単な,

(a)y = ax^{-b}

ベイズ因子やWAICで比べてみるとか。

my.prior <- prior(normal(0,3), nlpar = a) + prior(normal(0,1), nlpar = b)
fit2 <- brm(bf(m ~ (a*s^-b), a + b ~ 1, nl = TRUE),
	data = dat,
	prior = my.prior,
	save_all_pars=T)
summary(fit2)

f:id:kusanagik:20180914172915p:plain

f:id:kusanagik:20180914172933p:plain

bayes_factor(fit,fit2)
waic(fit, fit2)

fitがfit2より高かったけど…微妙。