草薙の研究ログ

英語の先生をやってます。

Rで繰り返し数列を振る

1, 2, 3, ... nという連番を振るなら，

1:10
seq(1,10,1)

などとやればよいのだけど，1からnまでの連番をm回繰り返したいときがある。
たとえば1, 2, 3, 1, 2, 3, 1, 2, 3みたいな。
これもrep関数で便利にできる。

rep(1:3, 3)

しかし，1, 1, 1, 2, 2, 2, 3, 3, 3みたいにしたいときめんどい。
for関数を使えばよい。でも，すぐfor関数を使ったりそれを埋め込んだりする男はディナーの支払いがスマートでない男のようなものだって祖母がいってた。
実はこれ，ちょっと気が効いたやり方がある。

#まずはさっきと同じ
rep(1:3, 3)
[1] 1 2 3 1 2 3 1 2 3

#これを行列にする
matrix(rep(1:3, 3),3,3)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3

#これを転置しちゃう
t(matrix(rep(1:3, 3),3,3))
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    3
[3,]    1    2    3

#これを戻す
as.numeric(t(matrix(rep(1:3, 3),3,3)))
[1] 1 1 1 2 2 2 3 3 3

うぇーい。できた。
などとこれまでずっと数年間この方法でやっていたら，最近になってrep関数にオプションがあることを知った。
「なんですかそのジジ臭いやりかたは！」と怒られた。ばあちゃん，ごめん。

#eachを使う
rep(1:3,each=3)
[1] 1 1 1 2 2 2 3 3 3

まぢで。
やはりなんでもすぐ調べる癖をもたなければならぬ。メロスは政治がわからぬらしいが私はなかなかスマートにはなれぬ。

全然わからない研究の話に巻き込まれたときに使う表現集

雑感 ライフハック

よくわからない研究の話に巻き込まれる！

私は，大学において外国語科目の授業を行い，外国語の教育，特に高等教育におけるそれに関することを一応研究上の専門としているのだけど，この分野は「〇〇学」と一括りにできるようなものではない。たとえば，私の学位は学術であって〇〇学ではないし，関連学会にはさまざまな「学」を名乗る方が混在している。まあ，平たくいえば，学際分野のひとつになっているわけ。

学際分野といえば聞こえはいいけど，要するにさまざまな経歴，背景，信条，技術をお持ちの方が近くにいっぱいいらっしゃるわけで，具体的にいえば，文学，言語学，教育学，心理学，最近では工学，経済学，社会学などに関連するキャリアをおもちの方々がひとつの学会に集まったり，査読をしあったりするというわけだ。上記のような名前が入る学会にこちらから勉強しにいくこともある。

根本的に私が不勉強なだけなのかもだけど，当然ながら，このような多様化が進むにつれて，コミュニケーションはうまくいかなくなってくる理屈だ。懇親会に勇気を振り絞って参加してみたものの，言葉を交わす相手がまったくなにをおっしゃっているかがわからないことが（少なくとも私には）ままある。はっきりいって専門用語が多すぎてついていけない。この分野にはアルファベット略の専門用語も多すぎる。なにか薬学とか化学の話をしているのかな？などと思ってしまう。

どうにかしてうまく乗り切りたい！

そのたびごとに，こちらが意を介していないような素振りを見せたり，露骨に不機嫌な顔をしたりするのは社会人として慎みたいところだ。…ああ，そうか，「カラオケに見られるような互恵主義的慣習がここでうまく機能しているのだなあ」とでも思いたい。私はあなたの話をまったく介していないし面白くもないし，お酒はもういいから帰ってホテルでhuluとかnetflixとか使って海外ドラマを鑑賞したいのだけど，ああ，それでも私の話を聞いてもらったのだから，今度はこちらも話を合わせるのが礼儀か，みたいな。つまり，そういうことみたい。

しかし，変に内容に質問したりコメントをしたりすると，学の浅いやつだと思われてしまう。できれば，まったくお話の内容には関知せずに，さもお話の内容を理解しているような態度を見せつつ，話し手を気持ちよくするような合いの手を入れたい。

合いの手を入れる！

そういうときに，以下のような表現が役に立つかも。5種類くらいをランダムにいっておけば間違いない。

  • …なるほど…そうなりますか…
  • それは先生らしい着眼点ですね
  • ほーほー，面白そうですね
  • なかなかそれをなさる研究者はいらっしゃらないでしょう
  • 本当に，〇〇を研究なさる方には頭が下がります
  • その研究テーマは，重要性に反して本当に研究が少ないですよね
  • まさに今の社会には必要な研究ですね
  • そのように思っている研究者は多いと思います
  • 測定の精度も問題というわけですね
  • それが何を測っているかを真剣に考えてみるわけですか
  • いろんなところにバイアスがあるのが科学ですからねえ
  • きっと素晴らしい論文になるでしょう
  • これまでのモデルでは明らかにうまくいかなかったですものね
  • そういう研究はこれまで聞いたことがありませんねえ（勉強不足で）
  • そういうところをしっかりしてないのが過去の研究の悪いところですね
  • それは今後流行ると確信していました
  • そのデータを取るのも簡単ではありませんよね
  • 私も実はずっとそうじゃないかと思っていました
  • 確かに問題がなかったとはいえませんものね

内容などに一切関知せずに，それっぽくなる。無情報がもっとも安全だというわけだ。コールドリーディングとかのテクニックに似ている。

話題を移行させる！

ちょっとだけ知識があるんだったら，不勉強をこちらから告白して御高説モードに移行させるのもひとつのテクだ。下手に質問するよりはるかに安全。

  • その分野は，〇〇（2010）の概論書を斜め読みしたぐらいで…
  • なかなか私のようなものにはその分野は手が出ないところにありまして
  • その分野は，難しくてなかなかついていけておりませんもので…

こうすると御高説モードになるから，いい塩梅で思い切り，

  • 先生，今日は勉強させていただいてありがとうございました。できれば今度，どこかでご講演をお願いしたいところですが，私などにはそんな裁量はなくて…えへへ

これでバッチリ。

ちょっと強引だけど，どうしても話を変えたいときは，

  • 先生，今回はどちらのホテルをお取りです？
  • せっかくですからなにか美味しいものを食べて帰りたいですねえ
  • 今年は晴れてよかったですねえ

などと，天気，土地，旅行など差し障りのないネタを入れて牽制する。中年以上のベテランの方には健康の話やご子息の話をしておけば安全。若手中堅は子育て，勤務，待遇，学会運営の話を，とにかく愚痴に共感すると小難しい研究の話をされずに済む。院生には業績や恋愛の話をしておく。このように話題転換も大事かな，と。

  • 先生はいつもご健康でいらっしゃいますねえ
  • ご子息はもう大学を出るころですか
  • お子さんはいくつになられるんでしたっけ？
  • 〇〇さんも大変そうですねえ
  • なかなか院生には大変な時期だよねえ
  • いきなり失礼だけど，えと，あの，聞きづらいんだけど…

これで話題が変えられるかも。これで学会で全然わからない研究の話に巻き込まれても安心だ。

 

…なんだこの記事？

 

私の一日とチャップマン＝コルモゴロフ方程式

r 統計

授業がない日のわたし，たとえば出張のない土日とか，10時から22時まで学校にいるとして，大抵は以下のような状態を自由に遷移している。

  • 研究関係のお仕事：論文の執筆，資料・論文・書籍を読む，データ分析，調査・実験の準備，やりたくない査読，共同研究者などへの連絡…
  • 事務手続きのお仕事：学内業務の書類，学会等の手続き，各種メール業務など…
  • 授業準備のお仕事：教材の準備，評価，学生支援など…
  • ネット徘徊：研究に関する情報収集，息抜きなど（SNSはやらない）
  • 喫煙休憩：ヘビースモーカーなので欠かせない，食事は抜くことが多い，昼寝はなかなかしない

ほとんどの場合，30分程度で気ままに別の状態にいったりきたりを繰り返している。これは定常性のあるマルコフ連鎖でモデル化できる。遷移確率表は，ざっとこんな感じ。

 

  研究 事務作業 授業関係 ネット 喫煙休憩

研究

 0.60 0.00 0.00 0.20 0.20
事務作業 0.00 0.30 0.30 0.00 0.40
授業関係 0.00 0.30 0.30 0.10 0.30
ネット 0.40 0.00 0.00 0.10 0.50
喫煙休憩 0.50 0.20 0.20 0.10 0.00

 

 研究に関することをやると，事務や授業関係にはいかず，繰り返しやり続けることが多く，ネットや喫煙をしないと頭が切り替わらない。ネットや喫煙がハブになっている。

問題は，もう少し長期的に見た場合，たとえば，1時間後，90分後，2時間後はどうなっているかということだ。今研究しているのだけど，1時間後，90分後，2時間後はなにをしているのか，または今喫煙しているのだけど，1時間後，90分後，2時間後はなにをしているのか。

これは有名なチャップマン＝コルモゴロフ方程式で予測できる。

 

チャップマン＝コルモゴロフ方程式 - Wikipedia

 

たとえば2ステップ後（1時間後）の遷移確率表はこうだ。

  研究 事務作業 授業関係 ネット 喫煙休憩
研究 0.54 0.04 0.04 0.16 0.22
事務作業 0.20 0.26 0.26 0.07 0.21
授業関係 0.19 0.24 0.24 0.07 0.26
ネット 0.53 0.10 0.10 0.14 0.13
喫煙休憩 0.34 0.12 0.12 0.13 0.29

 

3ステップ後（90分後）の遷移確率表はこうだ。

  研究 事務作業 授業関係 ネット徘徊 喫煙休憩
研究 0.50 0.07 0.07 0.15 0.22
事務作業 0.25 0.20 0.20 0.09 0.26
授業関係 0.27 0.20 0.20 0.10 0.24
ネット徘徊 0.44 0.09 0.09 0.14 0.25
喫煙休憩 0.40 0.13 0.13 0.12 0.22

 

4ステップ後（120分後）の遷移確率表はこうだ。

  研究 事務作業 授業関係 ネット 喫煙休憩
研究 0.47 0.08 0.08 0.14 0.22
事務作業 0.32 0.17 0.17 0.11 0.24
授業関係 0.32 0.17 0.17 0.11 0.24
ネット 0.44 0.10 0.10 0.14 0.22
喫煙休憩 0.40 0.12 0.12 0.13 0.23

 

時間が遅くなれば遅くなるほど，今何しているかにかかわらず，だんだん研究の方に向かっていくようだ（trivial...）。

 

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年もまた見てしまった…

はやい反応時間という表現。これが，外国語教育研究を含むある種の学術領域において，かなり慣習的な表現になっていることは知っている。たとえば，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などといえばいいかもだ。なんでも単位や形容詞が不安なときは，これが安パイ。

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

  • 反応が速い
  • 早く完了する

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

うん。

 

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

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

うっす。

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

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.