ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2017年12月13日

方向音痴の人は正しく「北」を指し示せるか──フォン・ミーゼス分布を使って角度データを分析する──

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の13日目のエントリー記事です。この記事では,角度や色相,あるいは1日の時間のような,循環する変数のモデリングに使えるフォン・ミーゼス分布を用いた分析について,実験心理学っぽい例を挙げて説明していきます。実際に分析するためのRとStanのコードも紹介します。

1. はじめに

 という訳で,アドカレデビューしてみました。改めまして,mutopsyこと武藤拓之です。大阪大学で認知心理学の研究をしています。いろんなものに興味がありますが,主にはヒトの空間認知過程について実験心理学の手法を用いて研究しています。空間認知の研究では変数として角度(30度とか270度とか)を用いることがよくあります。角度は循環する変数なので,360度回転すると一周して元の角度に戻ってきます。このような角度データに正規分布を当てはめてしまうと,循環する性質をうまく表現できず,例えば「30度付近のデータが得られる確率は390度付近のデータが得られる確率よりも高い」といった奇妙なことが起こってしまいます。そこでフォン・ミーゼス分布の出番という訳です。

2. 架空の実験データの例

 分布の説明をする前に,まずはデータの具体例を挙げたいと思います。以下の図1は架空の実験場面を描いたものです。実験室を訪れた実験参加者さんに,「北の方角を指差してください」と教示します。この実験室がある学部棟は独特な形をしており,実験参加者の正面(0度)から時計回りに120度の方向が北になっています。つまり,120度付近を正しく指差せたら「正解」となります。実験参加者が実際に指を差した方向(正面を0度とした時計回りの角度)をデータとして記録します。

vm_experiment.png
図1. 実験場面の例。実験室を訪れた実験参加者は,「北」を指差すように求められる。

 この実験を200人の実験参加者に対して行ったと想定します。そのうちの100人は方向感覚が良い人たちで,残りの100人は方向感覚が悪い人たち,すなわち方向音痴な人たちであったことが事前に分かっていたとします。このようにして得られた200人分のデータをプロットしたら,以下のような結果が得られました。なお,このデータはこの分析のために作った架空データです。このデータを生成するRコードは後述します。

vm_dataplot.png
図2. 指差し実験から得られた架空の角度データの分布。左右のグラフは同じデータを異なるグラフで表している。

 データを見てみると,方向感覚が良い人たちのグループは,多少の誤差はあれども,ほぼ正確に120度の方向を指し示せていることが読み取れます。一方,方向感覚が悪い人たちのグループの回答はかなりバラついています。バラついてはいるものの,均等に散らばっているというよりは,0度(=正面)付近を指す人の割合が若干高そうです。

3. フォン・ミーゼス分布について

 循環する変数が従う分布としては,フォン・ミーゼス分布がよく用いられます。フォン・ミーゼス分布は単峰かつ対称な分布で,循環正規分布とも呼ばれます。平均値パラメータμと集中度パラメータκ(カッパ) の2つのパラメータで形が決まります。集中度パラメータのκは非負の実数で,値が大きいほどμの付近に値が集中し,かつ正規分布に近づきます。κが0の時には一様分布となります。また,κは第1種変形ベッセル関数を用いることで,円周分散Vと円周標準偏差ν(ニュー)に変換できます。この円周標準偏差νは,正規分布の標準偏差σと同様の解釈ができます (豊田, 2017)。以下の図3は,μとκを変えた時にフォン・ミーゼス分布がどのような形状になるのかを示しています。ここでは0度と360度を両端とする確率密度関数を描画していますが,実際には横軸が循環していることに注意してグラフを見る必要があります。例えば中央上のグラフと右上のグラフは一見すると二峰分布に見えてしまうかもしれませんが,いずれも値がμの時にピークとなる単峰分布です。

vms.png
図3. いろいろなフォン・ミーゼス分布。

4. シミュレーション用のデータセットを作る

 以下に,フォン・ミーゼス分布から乱数を発生させて架空データを作るRコードを示します。2節で紹介した架空データはこのコードを用いて作成しました。循環する変数を扱うにはcircularパッケージが便利です。このパッケージのrvonmises関数を使うことでフォン・ミーゼス分布からの乱数を発生させることができます。また,ここでは詳しく説明しませんが,plot.circularという関数を使えば循環する変数を円周上にプロットすることができます。図2の左のグラフはこの関数を使って作成しました(ちなみに右のグラフは皆大好きggplot2で作成しました)。

Rコード01: フォン・ミーゼス分布に従う架空データを2セット作る
##前準備
library(circular)
rad_trns <- pi / 180 #これを度数に掛けるとラジアンに変換できる

##架空データのパラメータを設定
n1 <- 100    #群1のサンプルサイズ
n2 <- 100    #群2のサンプルサイズ

m1 <- 120 * rad_trns  #群1の平均値
k1 <- 50  #群1の集中度

m2 <- 0 * rad_trns  #群2の平均値
k2 <- 0.5  #群2の集中度

##円周分散と円周標準偏差を計算する
A1 <- besselI(k1, 1, expon.scaled=F)/besselI(k1, 0, expon.scaled=F)
v1 <- 1 - A1  #群1の円周分散
s1 <- sqrt(-2 * log(A1))  #群1の円周標準偏差

A2 <- besselI(k2, 1, expon.scaled=F)/besselI(k2, 0, expon.scaled=F)
v2 <- 1 - A2  #群2の円周分散
s2 <- sqrt(-2 * log(A2))  #群2の円周標準偏差

##乱数を生成
d1 <- as.numeric(rvonmises(n1,circular(m1),k1))  #群1のデータ
d2 <- as.numeric(rvonmises(n2,circular(m2),k2))  #群2のデータ
  • 前準備(1-3行目):circularパッケージとrstanパッケージを読み込み,rad_trnsという変数を作る。この変数rad_trnsを度数(e.g., 90度)に掛けるとラジアン(e.g., π/2)になり,この変数でラジアンを割ると度数になる。角度データを扱う時には作っておくと便利。ちなみにcircularパッケージには度数をラジアンに変換する関数radがあるのでこれを使っても良い。
  • 架空データのパラメータを設定(5-13行目):群1と群2のパラメータの真値とサンプルサイズを設定する。ここの値をいろいろ変えてみながら遊ぶと楽しい。
  • 円周分散と円周標準偏差を計算する(15-22行目):上で設定したパラメータを変形して,各群の真の円周分散と円周標準偏差を得る。besselI(x, n, expon.scaled=F) はn次の第1種変形ベッセル関数。
  • 乱数を生成(24-26行目):平均m1 or m2,集中度k1 or k2のフォン・ミーゼス分布に従う乱数をn1 or n2個生成する。rvonmises関数で引数のfromとtoを指定すれば,生成する乱数の範囲を指定できる (例えばfrom = -pi, to = pi; デフォルトではfrom = 0, to = 2 * pi)。また,デフォルトでは弧度法(ラジアン)のデータを渡すようになっているが,control.circular = list(units = "degrees")と指定すれば度数法のデータを渡すこともできる。詳しくはヘルプを参照。

5. Stanで分析を実行

 それでは,先ほどR上で生成した架空データにフォン・ミーゼス分布を当てはめて,パラメータを推定してみましょう。さらに,(1) 方向感覚が良い人たちと悪い人たちとの平均値差(最短角距離)はどのくらいか,(2) 方向感覚の悪い人たちの円周標準偏差は方向感覚が良い人たち円周標準偏差の何倍か,という2つの問いに答えるための生成量も推定してみましょう。以下にそのためのStanコードを示します。

Stanコード01: フォン・ミーゼス分布のパラメータを推定する
data {
  int N1;  //群1(高群)のサンプルサイズ
  int N2;  //群2(低群)のサンプルサイズ
  real Angle1[N1];  //群1の角度データ(ラジアン)
  real Angle2[N2];  //群2の角度データ(ラジアン)
}

parameters {
   real<lower=-pi(), upper=pi()> mu1;  //群1の平均
   real<lower=0, upper=2*pi()> mu2;  //群2の平均
   real<lower=0, upper=200> kappa[2];  //群1と群2のそれぞれの集中度
}

model {
  for (n1 in 1:N1){
    Angle1[n1] ~ von_mises(mu1, kappa[1]);  //群1のモデル
  }
  for (n2 in 1:N2){
    Angle2[n2] ~ von_mises(mu2, kappa[2]);  //群2のモデル
  }
}

generated quantities{
  real mu1_trns;  //mu1(ラジアン)を度数に変換
  real mu2_trns;  //mu2(ラジアン)を度数に変換
  real a[2];  //平均合成ベクトル長(円周分散・円周標準偏差の計算に使用)
  real v[2];  //群1と群2の円周分散
  real nu[2];  //群1と群2の円周標準偏差
  real nu_trns[2];  //nu(ラジアン)を度数に変換
  real mu_trns_diff;  //mu1_trnsとmu2_trnsの差(最短距離)
  real ratio_nu; //nu1とnu2の比

  #上で宣言した生成量を計算
  for (k in 1:2){
    a[k] = modified_bessel_first_kind(1,kappa[k])/modified_bessel_first_kind(0,kappa[k]);
    v[k] = 1 - a[k];
    nu[k] = sqrt(-2 * log(a[k]));
    nu_trns[k] = nu[k] / pi() * 180; 
  }
  mu1_trns = mu1 / pi() * 180;
  mu2_trns = mu2 / pi() * 180;
  mu_trns_diff = mu2_trns - mu1_trns;
  ratio_nu = nu[2] / nu[1];
}

  • dataブロック(1-6行目):Rから渡すデータを記します。
  • parametersブロック(8-12行目):推定するパラメータを宣言します。平均パラメータmu1とmu2の範囲については,レンジが2π (=360度) になるように上限と下限を指定しておきます。また,上限と下限の境界周辺にパラメータが多く存在する場合には収束効率が著しく低下し推定がうまくいかなくなります。これは,境界部分が計算上非連続な点になってしまうことに起因すると思われます。これを防ぐために,あらかじめデータの分布を見て平均値の目星を付け,その平均値がパラメータの範囲の中心に位置するように上限と下限を設定するのが良いでしょう。例えば,群2の平均値の真値は0度ですが,mu2の範囲を<lower = 0, upper = 2 * pi()>としてしまうと,パラメータ空間が「0度〜」の部分と「〜360度(= 2π)」の部分の二か所に分断されてしまい,0度(= 360度)に近い値のパラメータをうまくサンプリングできなくなってしまいます。また,この状態で得られたMCMCサンプルは左右両極にピークを持つ二峰分布となるため,そのまま平均値を計算すると,EAP(mu2) ≒ π(= 180度)となり,真値(0度)と乖離してしまいます(こちらに関しては,事後分布の要約としてMAP推定値と最高密度区間を用いることで対処できます)。こういった問題を避けるために,上記のコードではmu2の範囲を-pi()からpi()までに指定してあります。興味のある方は,この範囲をあえて<lower = 0, upper = 2 * pi()>に変えたらどうなるのかをぜひ試してみてください。トレースプロットやn_effを見ると上記の説明が腑に落ちると思います。
  • modelブロック(14-21行目):互いに独立な2つのフォン・ミーゼス分布から群1と群2のデータが生成されるモデルを記述します。Stanでは,von_mises(mu, kappa)でフォン・ミーゼス分布からの乱数を発生させることができます。引数のmuには度数ではなくラジアンを与えます。
  • generated quantitiesブロック(23-44行目):modified_bessel_first_kind(n,x)はn次の第1種変形ベッセル関数です。RのbesselI関数とは引数の順番が異なるので気を付けましょう。このブロックで生成しているxxx_trnsという変数は,ラジアンを度数に変換したものです。パラメータを解釈する時には度数の方が分かりやすいのでここで変換しています。mu_trns_diffはmu1_trnsとmu2_trnsの差の生成量ですが,角度は循環する変数なので,モデルを書く時にはパラメータの範囲や引き算の順序に気を配る必要があります。

 データをStanに渡すRコードは以下の通りです。分かりやすさのために,MCMCのパラメータを明示的に指定しています。

Rコード02: 角度データをStanに渡す
library(rstan)

##データをStanに渡す
data <- list(N1 = n1, N2 = n2, Angle1 = d1, Angle2 = d2)
stanmodel_vm <- stan_model(file="Stancode/von_mises.stan")
fit <- sampling(
  stanmodel_vm,
  data=data,
  seed=123,
  chains=4,
  iter=5500,
  warmup=500,
  thin=1
)

 このコードを走らせて得られたMCMCサンプルをsummary関数で要約すると以下のような結果が得られます。

von_mises_output1.png
図4. パラメータと生成量の推定結果をsummary関数で要約したもの。

 フォン・ミーゼス分布の平均パラメータ(今回の場合はmu1とmu2)の結果をsummary関数から読み取る際には注意が必要です。この出力結果に表示されているmean (=EAP) と分位点の値はMCMCサンプルから計算された値ですが,このMCMCサンプル自体は指定したlowerとupperの範囲におけるパラメータの集まりで,循環する性質を持っていません。例えば,真のμが0度付近の場合にパラメータの範囲を<lower = 0, upper = 2*pi()>と指定してしまうとEAP推定値は180度付近になってしまいますし,分位点に基づく95%ベイズ信頼区間に0度(360度)が含まれなくなります。そのため,フォン・ミーゼス分布の平均パラメータを要約する際には最初に事後分布のヒストグラムを目視し,それからMAP推定値と最高密度区間を計算するという手順を踏むのが良いでしょう。ヒストグラムのピークがグラフの両端に分断されてしまっている場合には,パラメータの範囲を指定し直して推定をやり直すのが賢明です。今回の分析では,μの真値がパラメータ範囲の中央付近となるように範囲をあらかじめ指定しておいたので,図4の結果をそのまま用いてもそれなりに真値と一致します。とはいえ基本的にはMAP推定値と最高密度区間を用いるほうが安全です。

 推定したパラメータと生成量の真値・MAP推定値・95%最高密度区間を以下の表1にまとめました。推定はうまくいっているようですが,群1に比べると群2の最高密度区間の幅が若干広めで,推定精度が相対的に低くなっています。これはもともとの真のκが小さいことに由来すると思われます。結果をまとめると次のようになります:(1) 方向感覚の良い人たちは北をほぼ正確に指差すことができ,個人差も比較的小さい。(2) 方向音痴な人たちには正面(0度)付近を北であると誤って回答する傾向がある。(3) 方向音痴な人たちの回答は個人差が大きく,円周標準偏差で見ると,方向感覚が良い人たちの8〜13倍ばらついている。(4) 方向感覚の良い人たちと悪い人たちの回答の平均値差(最短角距離)は100〜147度程度。※ただし,これらの結果はあくまでも架空データの結果なので,実際にこのような結論が得られるとは限りません。

表1. パラメータと生成量の真値・MAP推定値・95%最高密度区間
vm_summary.png

6. 結び

 この記事では,フォン・ミーゼス分布に従うデータをRとStanで分析する方法について解説しました。今回は分析に架空データを用いましたが,将来的には実際に実験で集めたデータを使って分析をしてみたいものです。個人的には色相データ(実験参加者に色を選ばせる)にフォン・ミーゼス分布を適用したら楽しいことになる予感がします。フォン・ミーゼス分布の混合分布なんかも面白そうです。フォン・ミーゼス分布は実験心理学とは割に相性の良い分布だとは思うのですが,いざ実際に研究で使おうと思うとなかなか難しいものです。面白い適用例などがありましたらTwitterなどで教えて頂けると嬉しいです。あと,明日もStan Advent Calendar 2017の記事を投稿します。正答率データを分析する方法(モデル)をベイズで比較するという趣旨の記事です。そちらもご覧になってくださると嬉しいです!→実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──

引用・参考文献

posted by mutopsy at 00:00 | 統計