2018-01-30

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

r 統計

もうしつこいのはわかっているけど，また，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

2018-01-30

「反応時間がはやい」のカテゴリー錯誤

雑感

2018年もまた見てしまった…

はやい反応時間という表現。これが，外国語教育研究を含むある種の学術領域において，かなり慣習的な表現になっていることは知っている。たとえば，google scholarで以下のような表現を検索してみると，かなりの件数が表示される。興味ある方は試して欲しい。

早い反応時間
速い反応時間
遅い反応時間
早い応答時間
速い応答時間
遅い応答時間
fast response time
fast reation time
slow response time
slow reaction time
fast response latency
slow response latency

日本語に限ったことでなくて，英語だと数千件ずつヒットするような具合だ。このような表現が慣習的に使用されていることは紛れもない事実だとしても，このような表現に私は違和感を感じる。そして，非常に誤解を招きやすい危険な表現だと思う。ちょっとここでこれに関連することを整理したい。

速さ，道のり，時間

速さ（speed）は，単位時間あたりの物体の位置の差（距離）だ。小学校で誰もが学習するように，速さ=道のり÷時間。反応時間は，名前がまさにそれを示すように時間なのだから，時間に対して速いという形容詞は，厳密にいえばカテゴリー錯誤だ。時間は長短という性質をもち，速さは速いか，遅いかだ。まちがった形容詞を名詞につけている。

反応時間の単位は，大抵の場合，秒か，ミリ秒。記録を終えた時刻から記録を始めた時刻を引いたものが反応時間。一方，速さの単位は，分子にある道のり。ある一定の時間あたりに，どれだけ位置に差が生じたかということ。

たとえば，道のりをこなす課題の量だと考えるとわかりやすい。単語が擬似語かを判断する課題をさせるとして，1分間あたりに50試行，これはまさに単位時間あたりのこなした量なので速さ。1分間あたりに50試行は，1分間あたりに25試行より速い。

一つ目のややこしさ：それぞれの関係

ただし，厄介なのが，道のり（に例えられる課題量や認知メカニズム？）が等しいと考えられるとき，そういう特殊な条件の下では，その道のりを移動し終える（課題を完遂する，認知メカニズムを終了する）までの時間が短ければ必然的に速く，その時間が長ければ必然的に遅い，といえるかもしれない。

日常言語であまり使い分けられているとはいえないけれど，速さ（speed）と速度（velocity）と早さ（earliness）は厳密にいえば異なることだと思うひともいる。学校でならったことによれば，速さはスカラーで，速度はベクトルだ。この区別はまだしも，早さは，ある事象の生じた時刻などが相対的に前のほうであるということだ。スタート時刻が一緒で，1kmを複数人が走るとしたとき，ほかのひとよりも早く着くひとがいる。このとき，その早く着いたひとは速く走るのかもしれないし，この条件下で着くまでの時間が短ければ早く着いて速く走るのかもしれない。非常にややこしい。落ち着いて考えたいところだ。しかし，早い反応時間というようないい方の違和感は変わりない。過去に自分もそのように書いてしまった論文があるのだけども。

二つ目のややこしさ：潜在変数と観測変数

もっともやっかいなことは，潜在変数と観測変数の関係，または操作化についてのこと。一般に，潜在変数としてなにかの技能や知識があって，反応時間はその潜在変数の影響を受ける観測変数として考えられることがある。たとえば，英語の読解力というものが少なくても心理測定的な含意のもとで存在するとされ，読解力（という名の潜在変数の値）が高ければ英語を読むのが速い（または読解力が高いひとは，そうでないひとよりも英語を読み終えるのが早い，または読解力が高いひとが示す読解時間は，そうでない人が示すそれよりも短い）といった関係性。または，かなり乱暴に操作化という手続きをして，反応時間は読解能力であると読み替えるかもしれない。この関係性や操作化による読み替えがカテゴリー錯誤の原因かもしれない。

たとえば，研究者のなかに，読解力＝反応時間というような置き換え規則があって，さらに，読解力が処理の速さ（processing speed）という側面をもつなどという理論的な先行研究を参照してたりしたら，読解力，つまり処理の速さは反応時間だから，反応時間が速い，というようないい方をついついしてしまうのかもしれない。つまり，読解力やそれに類するカテゴリーの性質を，反応時間という，まったく異なる性質をもつものに転移させているわけ。つまり，潜在変数がもつ性質と観測変数がもつ性質をごっちゃにしている。

できればこのようないい方を避けたい。これは論文全体の解釈にも影響を及ぼすことなので。

安全ないい方

既存の確立された学術領域における標準的な言葉遣いを避けることがよいことかはわからないけども（たぶんよくない），それでも，あくまでも反応時間に関しては，長短と書く方が無難かもしれない。

たとえば，（有意に長いとか短いとかはoutdatedだと思うけども）

条件Aにおける反応時間の平均値は，条件Bにおけるそれよりも有意に短い。
条件Aにおける反応時間の平均値は，条件Bにおけるそれよりも有意に長い。

または，もっとも安全で，そして汎用的な表現は，

条件Aにおける反応時間の平均値は，条件Bにおけるそれよりも大きな値を取った。
条件Aにおける反応時間の平均値は，条件Bにおけるそれよりも小さな値を取った。

これはこれで値が大小なのか高低なのかという問題があるけども，確率変数では最小値と最大値などともいうので，少なくとも変数の期待値を大小でいってもいいだろうと思う。英語では，greaterなどといえばいいかもだ。なんでも単位や形容詞が不安なときは，これが安パイ。

逆に，速いまたは早いを使うのだとしたら，

反応が速い
早く完了する

などといったいいかたの方が誤解を招かないかも。

うん。

2018-01-09

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

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

うっす。

2018-01-04

言語法則とベイズ的なカーブフィッティング

r 統計情報処理

続きもので書いていたつもりだったのだけど，放置しすぎて前回からの更新間隔がえげつないことになってしまった。前回は，

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

っていう話だった。えっと，「ある言語学的単位の平均長（y）は，その構成要素の長さ（x）の関数である」っていう古典的な観察があって，具体的には，この関数，

（a） $y = ax^{-b}e^{-cx}$

に従うだろう，という話だ。これをMezerath-Altmann法則と呼ぶんだそう（Altmann, 1980）。こういう言語法則系は話がとってもシンプルでいい。
たとえば，こんなデータ（Altmann, 1980）。

Morpheme Length（Syllables）	Syllable Length（Phonemes）
1	3.10
2	2.53
3	2.29
4	2.12
5	2.09

図にするとこんな感じ。確かに指数的減衰する関数ならフィットしそうだ。

f:id:kusanagik:20180104195751p:plain

普通に非線形最小二乗法でフィッティングしてもいいんだけど，世の中にはベイズ的なカーブフィッティングとかそういうのがあるそうで，RではFMEパッケージでそれができるんだそう。やってみよう。

これはメトロポリス法によって…とかそういう細かいところはぜんぶ置いといて（よくわからんし），まあこんな感じか。

#データ
s<-1:5
m<-c(3.1,2.53,2.29,2.12,2.09)

#モデルを書く
Model <- function(p, x) {
data.frame(x = x, 
N = p[1]*s^(-p[2])*exp(1)^(-p[3]*s))
}

Residuals <- function(p) {
(m - Model(p,s)$N)
}

#まずはパッケージの例通りにやってみる
P <- modFit(f = Residuals, p = c(0,0,0))
sP <- summary(P)
Covar <- sP$cov.scaled * 2.4^2/2
s2prior <- sP$modVariance

MCMC <- modMCMC(f = Residuals, 
p = P$par,jump = Covar,
niter = 5100,
var0 = s2prior,
wvar0 = NULL,
updatecov = 100,
burninlength=100)

#トレース図とMCMCの結果を可視化
par(mfrow=c(2,3))
title<-c("a","b","c")
for(i in 1:3){
plot(MCMC$par[,i],type="l",
xlab="Iteration",
ylab="",
main=title[i])
}

for(i in 1:3){
plot(density(MCMC$par[,i]),
type="n",
main=title[i])

polygon(density(MCMC$par[,i]),
col="gray98")
}

f:id:kusanagik:20180104201326p:plain

ほうほう。で，こんな感じになるというわけ。

f:id:kusanagik:20180104203102p:plain

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

2017-09-27

Menzerath-Altmann法則と一般化ガンマ分布（Eroglu, 2013）

統計

とても興味深いこと。Eroglu（2013）は，Menzerath-Altmann法則と一般化ガンマ分布の関係について説明している。

Eroglu, S. (2013). Menzerath-Altmann law for distinct word distribution analysis in a large text. Physica A, 392, 2775 – 2780

Menzerath-Altmann法則は，

（a） $y(x|A, b, c) = A x^b e^{-cx}$

という風に書ける。

ここで，b = 0，c ≒ 0のとき，

（b） $y(x|A,c) = Ae^{-cx}$

となる。これは普通の指数関数。

で，b ≒ 0，c = 0のとき，

（c） $y(x|A,b) = Ax^{-b}$

これはZipfの法則。

ここで，一般化ガンマ分布の確率密度関数を(d)のように書くとする（parameterizationもいろいろあるが…）。

（d） $y(x|α,β,θ) = \frac{β}{θγ(α)}(\frac{x}{θ})^{αβ-1}e^{-(\frac{x}{θ})^β}$

β = 1のとき，（d）式はガンマ分布の確率密度関数になるし，β=θのとき，ワイブル分布の確率密度関数になる。その上，ガンマ分布は，母数の値によっては，指数分布にもなるし，アーラン分布にもなる。

で，いちばん大事なんだけど， $A = (θ^αγ(α))^{-1}，b = α-1，c = \frac{1}{θ}$ のとき，（d）式は（a）式になるらしい（Eroglu, 2013, p. 2776）。なるほど。

…なんだか，とっても安心する。

2017-09-23

Menzerath-Altmann法則をRのnls関数で

r 統計教育

Menzerath-Altmann法則（またはMenzerath法則; e.g., Altmann, 1980; 基本的アイデアはMenzerathが発見，Altmannが定式化）とは，ある言語学的単位の増加が，その単位の構成要素の減少に帰着すること。またはその逆。

もっと一般的にいうと，Altmann（1980）は，「ある言語学的単位の長さは，その構成要素の長さの関数である」としている。たとえば，ある形態素長（形態素内の音節数などとして測定する）の値が高ければ，その形態素内の平均的な音節長（音節内の音素の数などとして測定される）の値は低くなるということ。ある言語において，1音節からなる形態素の平均的な音素数は3，2音節からなる形態素（ひとつつあたりの）の平均的な音素数は2.5，3音節からなる形態素の平均的な音素数は2…というように。

それで，基本的にはさまざまな単位，たとえば時間などでもこの法則が当てはまるとされ，実証研究も少なくないし，ゲノムとかタンパク質とか，そういう分野にも応用されているんだそうだ。で，その関係は，簡単な方程式で表されると。

もっとも一般的な式は，

（a） $y = ax^{-b}e^{-cx}$

ここでのyは構成要素（小さい方）の長さ，xは対象とする言語学的単位の長さ，a，b，cは自由母数。eはネイピア数。この式を使い，観測から，a，b，cを推定する，というわけだ。

一番単純な方法だと，最小二乗法を使用することができる。

Rには，非線形最小二乗法用の関数nlsがあるので，これを使えばいい。Altmann（1980）にインドネシア語の形態素と音節の関係についてのデータが示されている。これは，こんな感じ。

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

Morpheme Length（Syllables）	Syllable Length（Phonemes）
1	3.10
2	2.53
3	2.29
4	2.12
5	2.09

やってみる。

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

#可視化
plot(s,m,xlab="Morpheme Length",ylab="Syllable Length",pch=20,cex=2)

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

#サマリーなど
summary(fit)
AIC(fit)
BIC(fit)

#決定係数
cor(m,predict(fit))^2

#予測値をプロット
plot(s,m,xlab="Morpheme Length",ylab="Syllable Length",pch=20,cex=2)
points(s,predict(fit),col=2,pch=20)

f:id:kusanagik:20170923173803p:plain

赤がモデル理論値。もちろん原著の推定結果ともほぼ一致している。
このデータでは，以下のようになった。

（b） $y = 2.96x^{-0.36}e^{0.05x}$

この法則がどれだけ法則と呼ぶに相応しい一般化可能性をもつかといった面と，具体的に（a）式がどれだけ観測にフィットするかということは置いといて，ある単位と，その構成素の関係に一般的な法則をもとめる姿勢というのは，とても正当な考え方だと思う。

2017-09-20

ジニ係数とワイブル分布・対数正規分布・ガンマ分布の母数の関係

どうでもいいことなんだけど，ジニ係数はいくつかの分布であれば，分布の母数から直接的にもとまるんだそうだ。

ワイブル分布，対数正規分布，ガンマ分布を例に取ると，ワイブル分布とガンマ分布なら形状母数k（またはα）のみから，対数正規分布ならσのみから直接ジニ係数がもとまる。これらを計算するRの関数を作ってみる。

#ガンマ分布
gamma.gini<-function(a){
	gini<-gamma((2*a+1)/2)/k/gamma(a)/sqrt(pi)
	gini
}

#ワイブル分布
weibull.gini<-function(k){
	gini<-1-2^(-1/k)
	gini
}

#対数正規分布
lognorm.gini<-function(s){
	gini<-2*pnorm(s/2*sqrt(2))-1
	gini
}

あるデータの分布がワイブル分布に従っているとして，その形状母数kが2だったとしたら，

weibull.gini(2)

0.2928932

この関数にはベクトルを入れられるので，こんなグラフも作れる。

f:id:kusanagik:20170920145430p:plain

x<-seq(0.1,5,.01)
plot(x,weibull.gini(x),xlab="Parameter",ylab="Gini Coefficient",type="l")

ワイブル分布では，形状母数が大きくなればなるほどジニ係数は小さくなる。

…どうでもいいことなんだけど，ワイブル分布の形状母数について，MCMCのサンプルを得て，ジニ係数の信用区間を得ることもできるんじゃないだろうか。たとえば，

set.seed(1)
d<-rweibull(1000,3,1)
library(MCMCpack)

llf<-function(beta,x){
sum(log(dweibull(x,beta[1],beta[2])))
}

m<-MCMCmetrop1R(llf,theta.init=c(3,1),x=d)
gini<-weibull.gini(m[,1])
plot(density(gini),type="n",xlab="Gini",ylab="Density",main="")
polygon(density(gini),col=rgb(0,0,1,.2))

#信用区間
quantile(gini,c(.025,.975))


     2.5%     97.5% 
0.1963506 0.2140616

こんな感じで。

f:id:kusanagik:20170920150851p:plain

みたいな。