背景

いま一部の心理学領域で、認知モデルを立ててそのパラメタを推定する、あるいは、複数のモデルの良さを比較する、というデータ分析がスタンダードになりつつあります。また、その分析手法として、主にStanを使ったベイズ推定が取り上げられるようになっています。

しかし、一体どれだけの人（特に学部生）がベイズ推定についていけているでしょうか？　実際には、「ベイズ推定どころか最尤推定って何？　具体的にはどう計算するの？　どういうコードを書けばいいの？　そもそも認知モデルって何？」となってしまう¹場合が多いのでは？、というのが私の肌感覚です。

そういう「置いてけぼり感」を覚えてしまうことに無理はないです。現状、心理学（特に社会心理学）の授業で、認知モデルに関するトレーニングを受ける機会は少ないからです²。

このページの概要

そこでこのページでは、遅延割引課題を例として、認知モデルのパラメタをどうやって最尤推定するかを、具体的に見ていくことにします。推定にはRを使います。

なお、この課題については、『社会科学のためのベイズ統計モデリング』の第9章を参考にしています。より詳しく知りたい方は、そちらを読んでください。

遅延割引課題とは

遅延割引課題とは、「いますぐに小さい金額をもらうか、だいぶ先になっちゃうけどより大きい金額をもらうか」のどちらかを選ぶという課題です。

ここでは、「いますぐにr円をもらうか、n日後に50000円をもらうか」のどちらかを選ぶことになったとします。もっと具体的に考えるために、以下の2つの例を考えてみましょう。

まず、「いますぐに5000円をもらうか、30日後に50000円をもらうか」と聞かれたら、どっちを選びますか？　おそらく、多くの人は「30日後に50000円」を選ぶのではないでしょうか。

次に、「いますぐに45000円をもらうか、720日後に50000円をもらうか」と聞かれたら、どっちを選びますか？　ここでは逆に、「いますぐに45000円」を選ぶ人が多いと思います。

イメージはつきましたか？　この例からわかるように、遅延割引課題は「長期的利益のために、人々が目先の利益をどれくらい我慢できるか」を測定するものです。

遅延割引のモデル

さて、「人が長期的な利益をどれくらい良いと感じるか」は、以下の式（モデル）で表現できるとされています³：

\[ U(A,t) = U(A)\frac{1}{1+kt} \]

それぞれの記号は以下の内容を表しています：

• \(A\)：金額。ここでは50000円だと考えてください

• \(t\)：何日後に\(A\)（すなわち50000円）をもらえるか

• \(U(A, t)\)：\(t\)日後に\(A\)（50000円）をもらった場合の「嬉しさ」⁴

• \(U(A)\)：いますぐに\(A\)（50000円）をもらった場合の嬉しさ

• \(k\)：遅延割引パラメタ

では、この\(k\)（遅延割引パラメタ）とは具体的にどういうものなのでしょうか？　それを理解するために、\(k\)をいじってみましょう。

Figure 1の横軸は「t日」、縦軸は「そのとき50000円をもらったときの嬉しさ（すなわち、\(U(A,t)\)）」を表しています。ここでは、\(t\)日を0〜30で変化させ、\(k\)を0から0.5まで0.05刻みで変化させています。

この図からどういうことが読み取れるでしょうか？

library(tidyverse)

expand_grid(t = 0:30, k = seq(0, 0.5, 0.05)) %>%
  mutate(u = 50000 * (1 / (1 + k * t))) %>%
  ggplot(aes(x = t, y = u, group = k, color = k)) +
  geom_line() +
  scale_color_viridis_c() +
  labs(x = "t（日）", y = "t日後に50000円もらったときの嬉しさ") +
  theme_minimal(base_family = "ヒラギノ角ゴシック W3") +
  theme(axis.text = element_text(color = "#333333"))

まず、\(t = 0\)のとき、\(k\)がどんな値であれ、縦軸（50000円もらったときの嬉しさ）は50000円です。\(t = 0\)は「今すぐ」と同じだからです。

次に、\(t\)が大きくなるほど、「\(t\)日後に50000円」の嬉しさは低くなっています。

また、\(k\)が大きくなるほど（線が明るい色になるほど）、「\(t\)日後に50000円」の嬉しさは低くなっています。これを引っ張って考えると、「\(k\)が大きい人ほど、今すぐの利益に目がくらんでしまいがち」ということになります⁵。

2肢選択のモデル

では次に「2つの選択肢（目先の利益 vs. 長期的な利益）をどうやって選ぶか」のモデルを考えてみましょう。

ここでは「今すぐ5000円をもらうか、1ヶ月後に50000円をもらうか」を選ぶ状況を例とします。私たちは多くの場合、「1ヶ月後に50000円」を選ぶでしょう。

しかし、人間というものは、そこまで完璧ではありません。もし「ねえ、今すぐ5000円あげるけど、欲しいよね？」と100回も聞かれたら、何回かは誘惑にかられて「今すぐ5000円」を選んでしまうかもしれません。

このような現象を表すため、よくロジスティック関数⁶というものが用いられます。遅延割引課題の例では、以下のように書くことができます。

\[ P_{d} = \frac{1}{1+\exp(-\beta[U(A^{d})-U(A^{s})])} \]

記号の意味は以下のとおりです。

• \(P_{d}\)：遅延報酬（\(t\)日後に50000円）を選ぶ確率（0〜1の値）

• \(U(A^{d})\)：遅延報酬をもらった時の嬉しさ

• \(U(A^{s})\)：即時報酬（今すぐの報酬）をもらった時の嬉しさ

• \(\beta\)：逆温度パラメタ（嬉しさの差分にどれくらい敏感に反応するか）

詳しくは「ロジスティック関数」などで調べてほしいのですが、この式で大事なのは以下の2点です。

1. \(U(A^{d})-U(A^{s})\)が大きくなるほど\(P_{d}\)も大きくなる。すなわち、遅延報酬が即時報酬よりも魅力的であるほど、遅延報酬を選ぶ確率が高くなるということです。裏を返すと、即時報酬のほうが魅力的であるほど、遅延報酬を選ぶ確率が低くなる（つまり、即時報酬を選ぶ確率が高くなる）ということも表しています。

2. \(\beta\)が大きいほど、「嬉しさ」の差分に強く反応する。つまり、\(\beta\)が大きい場合、嬉しさの大きい方を忠実に選ぶようになるということです。逆に\(\beta\)が0に近いほど、どちらを選ぶかはランダムに近づく（「嬉しさ」の影響力が弱くなり、選択肢をでたらめに選ぶようになる）ということも表しています。

---

さて、遅延割引課題を用いた研究で知りたいのは、ある人（参加者）がどれくらい目先の利益に目がくらみがちか、そして、ある人（参加者）がどれくらい「嬉しさ」の差分に強く反応するかということです。これを言い換えると、ある人がどれくらいの\(k\)（遅延割引パラメタ）を持っているか、ある人がどれくらいの\(\beta\)（逆温度パラメタ）を持っているかということになります。

認知モデルのパラメタを推定するというのは、まさに「ある人がどれくらいのパラメタ（認知的な傾向）を持っているか」を探ろうとする分析を意味しています。

実験データを作ってみる

ここで具体的な分析方法を説明したいのですが、まだ手元にはデータがありません。

そこで、ある1人の参加者が50回の遅延割引課題に取り組んだと仮定します。具体的には「\(t\)日後に50000円をもらうか、今すぐ\(r\)円をもらうか」を選んだとしましょう。各回で、日数\(t\)は｛30, 90, 180, 360, 720｝のいずれか、今すぐの報酬\(r\)は｛5000, 1000, …, 50000｝のいずれかとします（5通り×10通りなので50回の課題、ということです）。では、下のコードのように仮想的なデータを作り、dataというデータフレームとして保存してみましょう⁷。

set.seed(1)
data <- expand_grid(
  t = c(30, 90, 180, 360, 720),
  r = seq(5000, 50000, 5000)
  ) %>%
  mutate(
    u_delay = 50000 * (1 / (1 + 0.01 * t)),
    p = 1 / (1 + exp(-0.00005 * (u_delay - r)))
  ) %>%
  group_by(row_number()) %>%
  mutate(choice = rbinom(1, 1, p))

では、この参加者は遅延報酬（「\(t\)日後に50000円」）を何回選んだのでしょうか？　遅延報酬を選んだ場合、dataのchoiceという変数は1になっています。逆に即時報酬（「すぐに\(r\)円」）を選んだ場合、choiceは0です。

参加者の選択を下にプロットしてみました。横軸は遅延日数\(t\)、縦軸は選択（0か1）、各パネルの上の数字は即時報酬の金額を表しています。

data %>%
  ggplot(aes(x = t, y = choice)) +
  geom_point() +
  scale_y_continuous(breaks = c(0, 1)) +
  facet_wrap(~r)

素朴にパラメタを推定してみる

さて、どのようにパラメタを推定すればよいのでしょうか？　まず、素朴なアイデアから出発してみます。

たとえば、「なんとなくだけど、\(k\)は0.1で、\(\beta\)は0.0001じゃね？」と考えたとしましょう。このとき、これらのパラメタの値はどれくらいふさわしいでしょうか？　また、そのふさわしさはどのように評価（計算）されるべきでしょうか？

とりあえず、上で説明した2つのモデルに代入してみましょう。ここでは、「今すぐ5000円」と「30日後に50000円」を比べた場合を考えてみます。

t <- 30
k <- 0.1
beta <- 0.0001

# u_delay：遅延報酬のU
u_delay <- 50000 * (1 / (1 + 0.1 * t))

# 30日後に50000円を選ぶ確率
p <- 1 / (1 + exp(-beta * (u_delay - 5000)))

print(p)

## [1] 0.6791787

さて、\(k = 0.1\)、\(\beta = 0.0001\)としたとき、遅延報酬を選ぶ確率は0.6791787になりました。では、実際のデータはどうだったでしょうか？ つまり、「今すぐ5000円 vs. 30日後に50000円」のとき、どっちを選んでいたでしょうか？

Figure 2を見ると、参加者が遅延報酬（30日後に50000円）を選んでいたということがわかります。では、上の計算で出てきた「遅延報酬を選ぶ確率は0.6791787」と「参加者が実際に遅延報酬を選んでいたという事実」はどのように関係しているのでしょうか。

「遅延報酬を選ぶ確率」というのは、「『遅延報酬を選ぶ』という結果が得られる確率」と言い換えることができます。この「『遅延報酬を選ぶ』という結果が得られる確率」は、パラメタ\(k\)と\(\beta\)を色々な値にすれば、それに応じて様々な値に変動します。

ここで、「参加者が実際に遅延報酬を選んで」いたとしましょう。このとき、適当にパラメタを代入して、「『遅延報酬を選ぶ』という結果が得られる確率」が高かった場合、それは何を意味しているでしょうか。

日常言語で表すと、「参加者は遅延報酬を選んでいた。適当にパラメタを代入したら、遅延報酬を選ぶ確率は高いらしい。これって実際のデータと結構近いじゃん」ということになります。さらにこれを引っ張ると、「結構近いんだから、このパラメタって、もしかして良い線行ってるんじゃね？」となります。逆に、予測が外れた場合は、そのパラメタはあまりよろしくないと言えます。

つまり、「ある結果が得られている。そこで、適当にパラメタを設定したら、その結果が得られる確率も高いようだ。だったら、そのパラメタは『ふさわしい』んじゃないか」ということです。このようなパラメタの「ふさわしさ」（パラメタを設定したとき、手元にあるデータが得られる確率）を真面目に言うと「尤度」になります。また、尤度をもとにパラメタを推定する方法を最尤推定と言います。

この例では、「いま、参加者が遅延報酬を選んだことがわかっている。そこで、\(k = 0.1\)、\(\beta = 0.0001\)としたとき、遅延報酬を選ぶ確率は0.6791787だった」ということになります。すなわち、尤度は0.6791787です⁸。

尤度を計算する

では、\(k = 0.1\)、\(\beta = 0.0001\)とし、全50回の尤度を計算してみましょう。まず、以下のコードを実行します。

t <- 30
k <- 0.1
beta <- 0.0001

data %>%
  # 遅延報酬のUと、遅延報酬を選ぶ確率pを計算
  mutate(
    u_delay = 50000 * (1 / (1 + k * t)),
    p = 1 / (1 + exp(-beta * (u_delay - r)))
  ) %>%
  # 遅延報酬を選ぶ確率pを抽出
  pull(p)

##  [1] 0.679178699 0.562176501 0.437823499 0.320821301 0.222700139 0.148047198
##  [7] 0.095349465 0.060086650 0.037326887 0.022977370 0.500000000 0.377540669
## [13] 0.268941421 0.182425524 0.119202922 0.075858180 0.047425873 0.029312231
## [19] 0.017986210 0.010986943 0.441064710 0.323695074 0.224986141 0.149714490
## [25] 0.096490497 0.060834074 0.037802587 0.023274618 0.014247244 0.008690106
## [31] 0.409782432 0.296323939 0.203450772 0.134137019 0.085891464 0.053918000
## [37] 0.033411753 0.020535219 0.012556698 0.007653837 0.393766567 0.282619107
## [43] 0.192864008 0.126583888 0.080801479 0.050617863 0.031325176 0.019236783
## [49] 0.011756686 0.007163930

上に表示されているのは、遅延報酬を選ぶ確率\(p\)です。

しかし、これをそのまま最尤推定に使うことはできません。なぜなら、尤度とは「あるパラメタを設定したときに、手元にあるデータが得られる確率」だからです。実験では、参加者が即時報酬を選ぶ場合も当然あるので、その場合は「即時報酬を選ぶ確率」を使わないといけません。したがって、各回の尤度（likelihood）は以下のようになります。

likelihood <- data %>%
  # 遅延報酬のUと、遅延報酬を選ぶ確率pを計算
  # さらに、参加者の選択に応じて条件分岐をして尤度を計算
  mutate(
    u_delay = 50000 * (1 / (1 + k * t)),
    p = 1 / (1 + exp(-beta * (u_delay - r))),
    likelihood = if_else(choice == 1, p, 1 - p)
  ) %>%
  pull(likelihood)

likelihood

##  [1] 0.67917870 0.56217650 0.43782350 0.67917870 0.22270014 0.85195280
##  [7] 0.90465054 0.06008665 0.03732689 0.97702263 0.50000000 0.37754067
## [13] 0.73105858 0.18242552 0.88079708 0.92414182 0.04742587 0.02931223
## [19] 0.98201379 0.01098694 0.55893529 0.32369507 0.77501386 0.85028551
## [25] 0.90350950 0.93916593 0.96219741 0.97672538 0.01424724 0.99130989
## [31] 0.40978243 0.70367606 0.79654923 0.86586298 0.08589146 0.94608200
## [37] 0.03341175 0.97946478 0.98744330 0.99234616 0.60623343 0.28261911
## [43] 0.19286401 0.87341611 0.91919852 0.05061786 0.96867482 0.98076322
## [49] 0.98824331 0.99283607

対数尤度を計算する

では、ここで全50回の尤度を1つの指標にまとめたいと思います（まとめたいですよね？）。このとき、直観的には尤度を足し算すれば良いんじゃないか、と思うかもしれませんが、それは間違いです。

正解は「かけ算」です。詳しい説明は適宜教科書を読むなり、ググってもらえるとありがたいです⁹。簡単な例で言えば、「3回連続でじゃんけんに勝つ確率は？」と聞かれたとき、「1/3を3回足して1」ではなく、「1/3を3回かけ算して1/27」と答えるのが正解になるのと同じ理屈です。

では、上で求めた尤度を全部かけ算してみましょう。

# prod()で、要素を全部かけあわせることになる
prod(likelihood)

## [1] 1.143327e-20

これは0.00000000000000000001143327を表しています。

ちっさ！って思うのではないでしょうか。と同時に、「やっぱ、パソコンはこんな細かい計算できて偉いなあ」と思う人がいるかもしれません。

しかし、残念ながらそれは間違いです。あまりに数字が小さいと、パソコンの計算にも誤差が生じてしまいます¹⁰。

これを防ぐため、尤度を計算するときは、対数を取るというのが決まりになっています。なぜなら、かけ算の対数を取ると足し算になるという便利な性質があるからです¹¹。対数を取って求めた尤度を「対数尤度」と言います。

また、対数は単調増加関数なので、「尤度の積が大きい（もっともらしい）場合は、その対数も大きくなる」という関係が成り立ちます。対数、便利です。

では、対数尤度を計算し、全部を足し算してみましょう。

sum(log(likelihood))

## [1] -45.91776

確認のため、尤度を全部掛け算したやつの対数を取ってみましょう。

log(prod(likelihood))

## [1] -45.91776

ちゃんと2つが一致しました。ただし、データの数が多いほどかけ算の誤差は大きくなるので、やはり対数を取ってから足し算をするようにしましょう。

optimでパラメタを推定する

さて、\(k = 0.1\)、\(\beta = 0.0001\)としたときの50回の対数尤度は-45.91776になりました。

で？　それが何？

そうです。これだけでは、まだ何も言えていません。他のパラメタの候補も調べないと、どのパラメタが良さそうかはわかりません。

では、どうすればよいでしょうか？　素朴に考えたら、有り得そうな\(k\)や\(\beta\)を片っ端から代入して対数尤度を計算する、というのが良さそうです。しかし、そんな時間は私たちにはありません。

そこで使うのが、Rの最適化関数optimです。最適化関数とは、ある関数の最大値／最小値を求めるための道具だと考えてください。

いま私たちがやりたいのは、対数尤度が最も大きくなりそうなパラメタ\(k\)と\(\beta\)の組み合わせを調べるということです。つまり、optimを使って、対数尤度（関数）が最大化されるような\(k\)と\(\beta\)を見つければ目的達成ということになります。

では、具体的な書き方を説明します。まず、対数尤度関数（ll_delay）を自分で定義します。

# 関数を定義する
ll_delay <- function(param, data) {

  # param[1]はk、param[2]はbetaを表している
  # dataは参加者のデータフレームそのものを表している

  data %>%
    # 毎回の対数尤度を計算する
    mutate(
      u_delay = 50000 * (1 / (1 + param[1] * t)),
      p = 1 / (1 + exp(-param[2] * (u_delay - r))),
      likelihood = if_else(choice == 1, p, 1 - p),
      ll = log(likelihood)
    ) %>%
    # 最後に全部足し合わせる
    pull(ll) %>%
    sum()
}

これをoptimに突っ込みます。parには推定したいパラメタの初期値、fnには最適化したい関数を入れます。control = list(fnscale = -1)とすると、関数の最大化が行われるようになります。その他の引数（この例ではdata）では、ll_delayにどんな引数を渡すかを指示してあげます。

optim(par = c(0.1, 0.0001), fn = ll_delay, data = data, control = list(fnscale = -1))

## $par
## [1] 0.0132120062 0.0000423154
## 
## $value
## [1] -30.70504
## 
## $counts
## function gradient 
##       89       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

結果のparが推定値を表しています。\(k\)が0.0132120062、\(\beta\)は0.0000423154でした。

では、正解は何だったのでしょうか？　正解は、データを作ったコードの中に隠れています。改めてコードを見てみましょう。

set.seed(1)
data <- expand_grid(
  t = c(30, 90, 180, 360, 720),
  r = seq(5000, 50000, 5000)
  ) %>%
  mutate(
    u_delay = 50000 * (1 / (1 + 0.01 * t)),
    p = 1 / (1 + exp(-0.00005 * (u_delay - r)))
  ) %>%
  group_by(row_number()) %>%
  mutate(choice = rbinom(1, 1, p))

最初は詳しく説明しませんでしたが、このコードでは、遅延割引と2肢選択のモデルを用いて、参加者の選択を仮想的に作っていました。とくに7行目と8行目を見てください。ここでは、参加者が\(k=0.01\)、\(\beta=0.00005\)というパラメタを持つように設定しています。さて、推定結果は\(k=0.0132120062\)、\(\beta=0.0000423154\)でした。まあまあ良い線行ってるのではないでしょうか。