スタートページ＞ JavaScript＞ 他言語＞ Ｒ 目次

ベイズ推論

ベイズ推定と単純な例題

ベイズ推定（Bayesian inference）とは、事前に設定された確率データ（事前確率行列）があり、そのうちの幾つかが実現したとき、そをの起因事象を推定する（事後確率を求める）という確率の方法で、ベイズの定理に基づいています。

第１表のような事前確率行列があります。例えば、都会ではネコが首輪をしている確率は 0.7、イヌを木の上で見かける確率は 0.1、サルが二足直立をしている確率は 0.8 であることが、何らかの調査でわかっている（どんな調査？）とします。

あるとき、首輪をしている動物を木の上にいるのを見かけました（これが実現特性です）。この動物が何であるかを統計的に知りたい（事後確率）のです。

実現特性が一つのとき

「首輪をしていた」という事実だけがあったとします。第１表の「首輪」の列だけが対象になります。これだけでも最もありそうなのはイヌだとわかります。
しかし、事後確率という考え方では「首輪をしていた」とは、この列の合計行の値が「確率＝１」になったのですから、そうなるように各列を合計で割る操作をします（第２表）。

実現特性が複数のとき

「首輪をしている動物が木の上にいるのを見た」場合です。
「首輪～」と「木の上～」は互いに独立の事象ですから、「首輪～、かつ、木の上～」の確率は、両者の積になります。
そして、「実現特性が一つのとき」と同様に合計を１になるように修正すると第３表になります。
これより、
　　ネコ　５６.５％
　　サル　２９.０％
などが得られます。

Ｒによるプログラム

# ============ 入力データ
事前確率行列 <- matrix(c(0.7, 0.5, 0.2,   # ネコ
                         0.8, 0.1, 0.3,   # イヌ
                         0.2, 0.9, 0.8,   # サル
                         0.1, 0.1, 0.6),  # クマ
                   nrow=4, byrow=T)
実現特性 <- c(1, 2)                       # 1（首輪）、2（木の上）を見た

# ============= 事後確率の計算
# ====== 作業変数
行数 <- nrow(事前確率行列)
列数 <- ncol(事前確率行列)

合計行 <- matrix(0, nrow=1, ncol=列数)
事後確率列 <- matrix(1, nrow=行数, ncol=1)

# ====== 計算
for (特性列 in 実現特性){
    合計行[特性列] <- sum(事前確率行列[ , 特性列])
    事後確率列 <- 事後確率列 * 事前確率行列[ , 特性列] / 合計行[特性列]
}
事後確率列合計 <- sum(事後確率列)
事後確率列 <- 事後確率列 / 事後確率列合計
事後確率列

# ============== 結果出力
順位 <- order(-事後確率列)
sprintf("%d, %8.3f", 順位[1],事後確率列[順位[1]])
sprintf("%d, %8.3f", 順位[2],事後確率列[順位[2]])

# 事後確率列
    # 0.56451613
    # 0.12903226
    # 0.29032258
    # 0.01612903
# 結果出力
    1   0.565'
    3   0.290'

「事後確率の計算」の関数化

上のプログラムで「事後確率の計算」がこのロジックの中核であり、入力データの影響を受けないので、関数にしてみました。この事後確率計算関数をライブラリにしておけば、この部分を記述する必要はありません。

# ============= 事後確率の計算
事後確率計算関数 <- function(事前確率行列, 実現特性) {
    行数 <- nrow(事前確率行列)
    列数 <- ncol(事前確率行列)
    合計行 <- matrix(0, nrow=1, ncol=列数)
    事後確率列 <- matrix(1, nrow=行数, ncol=1)

    for (特性列 in 実現特性){
        合計行[特性列] <- sum(事前確率行列[ , 特性列])
        事後確率列 <- 事後確率列 * 事前確率行列[ , 特性列] / 合計行[特性列]
    }
    事後確率列合計 <- sum(事後確率列)
    事後確率列 <- 事後確率列 / 事後確率列合計
    return(事後確率列)
}

# ============ 入力データ　（これ以降を記述するだけでよい）
事前確率行列 <- matrix(c(0.7, 0.5, 0.2,
                         0.8, 0.1, 0.3,
                         0.2, 0.9, 0.8,
                         0.1, 0.1, 0.6),
                   nrow=4, byrow=T)

# ============ 関数の利用
事後確率列 <- 事後確率計算関数(事前確率行列,  c(1,2))　
事後確率列

順位 <- order(-事後確率列)
sprintf("%d, %8.3f", 順位[1],事後確率列[順位[1]])
sprintf("%d, %8.3f", 順位[2],事後確率列[順位[2]])

実務適用での問題点

このモデルは単純なわりに、適用できる分野が広く、実用性が高いように思われます。
ところがこの 事後確率計算関数 を実務で利用するのは不適切です。

特性間の独立性への考慮

このモデルでは「各特性が独立であると仮定」しますが、「木の上にいるときは二足直立になりやすい」というように、特性間に相関があることが多いのです。それを排除するには共分散分析などが必要です。実務では、むしろ「Ａ特性が発生したときのＢ特性の発生確率」のようなデータが重要なこともあります。

計算精度への考慮

このモデルでは多くの数の掛け算が行われます。事前確率は全て１より小さい値ですから、すぐにアンダーフローが生じ、誤差が大きくなります。それで log にして加算する手段が必要になります。すると合計を求めるのが掛け算結果になってしまうので、合計は元のデータでの加算のほうが適切です。このように、数値計算での工夫が必要になります。

これらの対策を講じたモデルが多数公開されていますが、その内部構造はかなり複雑なようです。