はじめに

　機械学習でも登場する情報理論におけるエントロピーを1つずつ確認していくシリーズです。

　この記事では、結合エントロピーを扱います。

【前の内容】

【今回の内容】

はじめに
結合エントロピーの定義
- 定義
- 可視化
参考文献
おわりに

結合エントロピーの定義

　結合事象(joint event)の平均情報量(avarage information)である結合エントロピー(joint entropy)の定義を数式とグラフで確認します。エントロピーについては「平均情報量(エントロピー)の定義 - からっぽのしょこ」、条件付きエントロピーについては「条件付きエントロピーの定義 - からっぽのしょこ」を参照してください。

定義

　まずは、結合確率(joint probability)と結合エントロピーの定義を数式で確認します。結合確率は同時確率(simultaneous probability)とも呼ばれます。

結合確率の定義式

　事象 $x, y$ に依存関係がないとき、結合確率(結合事象の生起確率)は、「各事象の確率」の積で定義されます。

$\displaystyle p(x, y) = p(x) p(y)$

　依存関係がある場合は、「条件の確率」と「条件付き確率」の積で定義されます。

$\displaystyle p(x, y) = p(x) p(y | x)$

　依存関係がある場合については「条件付きエントロピーの定義」で扱います。

結合エントロピーの定義式

　事象 $x$ の自己情報量 $I(x)$ は、生起確率 $p(x)$ の負の対数で定義されるのでした。定数の底は $2$ (またはネイピア数 $e$ )とします。

$\displaystyle I(x) = - \log_2 p(x)$

　2つの事象系 $\mathbf{x} = (x_1, \cdots, x_n), \mathbf{y} = (y_1, \cdots, y_m)$ の確率分布を $p(\mathbf{x}), p(\mathbf{y})$ とします。
　 $\mathbf{x}, \mathbf{y}$ が独立なとき、 $\mathbf{x}, \mathbf{y}$ の結合エントロピー $H(\mathbf{x}, \mathbf{y})$ は、「結合確率の自己情報量」の期待値で定義されます。

$\displaystyle \begin{align} H(\mathbf{x}, \mathbf{y}) &= \mathbb{E}_{p(\mathbf{x}, \mathbf{y})}[I(x_i, y_j)] \\ &= \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) I(x_i, y_j) \\ &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) \log_2 p(x_i, y_j) \tag{1} \end{align}$

　「結合確率」と「結合確率の自己情報量(負の対数)」の積和で計算できます。
　多変数のエントロピー(平均情報量)と言えます。

　生起確率 $p(x)$ と自己情報量 $I(x)$ が非負の値 $p(x) \geq 0, I(x) \geq 0$ なので、結合エントロピーも非負の値 $H(\mathbf{x}, \mathbf{y}) \geq 0$ をとります。エントロピーも非負の値 $H(\mathbf{x}) \geq 0$ です。

他のエントロピーとの関係

　続いて、 $\mathbf{x}, \mathbf{y}$ に依存関係がある場合とない場合でそれぞれ式(1)を整理して、エントロピー(平均情報量)と条件付きエントロピーとの関係を導出します。

依存関係がない場合

　 $x_i, y_j$ に依存関係がないとして結合確率を分割 $p(x_i, y_i) = p(x_i) p(y_j)$ します。

$\displaystyle \begin{align} H(\mathbf{x}, \mathbf{y}) &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) \log_2 p(x_i, y_j) \tag{1}\\ &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j) \log_2 \Bigl( p(x_i) p(y_j) \Bigr) \\ &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j) \Bigl( \log_2 p(x_i) + \log_2 p(y_j) \Bigr) \\ &= - \sum_{i=1}^n \sum_{j=1}^m \Bigl\{ p(x_i) p(y_j) \log_2 p(x_i) + p(x_i) p(y_j) \log_2 p(y_j) \Bigr\} \end{align}$

　対数の性質より $\log_a (x y) = \log_a x + \log_a y$ です。
　さらに、総和についての括弧を展開します。

$\displaystyle \begin{aligned} H(\mathbf{x}, \mathbf{y}) &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j) \log_2 p(x_i) - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j) \log_2 p(y_j) \\ &= - \Biggl\{ \underbrace{ \sum_{j=1}^m p(y_j) }_{1} \Biggr\} \sum_{i=1}^n p(x_i) \log_2 p(x_i) - \Biggl\{ \underbrace{ \sum_{i=1}^n p(x_i) }_{1} \Biggr\} \sum_{j=1}^m p(y_j) \log_2 p(y_j) \\ &= - \sum_{i=1}^n p(x_i) \log_2 p(x_i) - \sum_{j=1}^m p(y_j) \log_2 p(y_j) \end{aligned}$

　 $p(x_i)$ は $j$ の影響を受けないので $\sum_j$ の外に出せます。また、 $p(x_i)$ の $i$ に関する総和は、全事象の確率の和なので、 $\sum_i p(x_i) = 1$ です。同様に、 $p(y_j)$ は $i$ の影響を受けないので $\sum_i$ の外に出せ、 $\sum_j p(y_j) = 1$ です。
　前の因子は $\mathbf{x}$ のエントロピー、後の因子は $\mathbf{y}$ のエントロピーの定義式なので、次の関係が得られます。

$\displaystyle H(\mathbf{x}, \mathbf{y}) = H(\mathbf{x}) + H(\mathbf{y}) \tag{2}$

　事象間で依存しないとき「結合エントロピー」は、各事象の「エントロピー」の和です。

　エントロピーは非負の値なので、式(2)から、結合エントロピーは各事象のエントロピー以上の値なのが分かります。

$\displaystyle \begin{aligned} H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{x}) \\ H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{y}) \end{aligned}$

依存関係がある場合

　 $x_i$ を条件として結合確率を分割 $p(x_i, y_i) = p(x_i) p(y_j | x_i)$ して、式を整理します。

$\displaystyle \begin{align} H(\mathbf{x}, \mathbf{y}) &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) \log_2 p(x_i, y_j) \tag{1}\\ &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j | x_i) \log_2 \Bigl( p(x_i) p(y_j | x_i) \Bigr) \\ &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j | x_i) \Bigl( \log_2 p(y_j | x_i) + \log_2 p(x_i) \Bigr) \\ &= - \sum_{i=1}^n \sum_{j=1}^m \Bigl\{ p(x_i) p(y_j | x_i) \log_2 p(y_j | x_i) + p(x_i) p(y_j | x_i) \log_2 p(x_i) \Bigr\} \end{align}$

　さらに、総和についての括弧を展開します。

$\displaystyle \begin{aligned} H(\mathbf{x}, \mathbf{y}) &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j | x_i) \log_2 p(y_j | x_i) - \sum_{i=1}^n \sum_{j=1}^m p(x_i) p(y_j | x_i) \log_2 p(x_i) \\ &= - \sum_{i=1}^n \Bigl\{ p(x_i) \sum_{j=1}^m p(y_j | x_i) \log_2 p(y_j | x_i) \Bigr\} - \Biggl\{ \underbrace{ \sum_{j=1}^m p(y_j | x_i) }_{1} \Biggr\} \sum_{i=1}^n p(x_i) \log_2 p(x_i) \\ &= - \sum_{i=1}^n \Bigl\{ p(x_i) \sum_{j=1}^m p(y_j | x_i) \log_2 p(y_j | x_i) \Bigr\} - \sum_{i=1}^n p(x_i) \log_2 p(x_i) \end{aligned}$

　 $p(y_j | x_i)$ は $i$ の影響を受けないので $\sum_i$ の外に出せます。 $p(y_j | x_i)$ の $j$ に関する総和は、全事象の確率の和なので、 $\sum_j p(y_j | x_i) = 1$ です。
　前の因子は(結合確率の形に戻すと)条件付きエントロピー、後の因子はエントロピーの定義式なので、次の関係が得られます。

$\displaystyle \begin{align} H(\mathbf{x}, \mathbf{y}) &= - \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) \log_2 p(y_j | x_i) - \sum_{i=1}^n p(x_i) \log_2 p(x_i) \\ &= H(\mathbf{y} | \mathbf{x}) + H(\mathbf{x}) \tag{3} \end{align}$

　事象間で依存するとき「結合エントロピー」は、「条件付きエントロピー」と条件の「エントロピー」の和です。

　また、条件付きエントロピーについて整理すると、次の関係が得られます。

$\displaystyle H(\mathbf{y} | \mathbf{x}) = H(\mathbf{x}, \mathbf{y}) - H(\mathbf{x})$

　「条件付きエントロピー」は、「結合エントロピー」と条件の「エントロピー」の差です。この関係は、条件付きエントロピーの定義式からも導出できます。

　条件を入れ替えた場合も同様に求められます。

$\displaystyle \begin{align} H(\mathbf{x}, \mathbf{y}) &= H(\mathbf{x} | \mathbf{y}) + H(\mathbf{y}) \tag{4}\\ H(\mathbf{x} | \mathbf{y}) &= H(\mathbf{x}, \mathbf{y}) - H(\mathbf{y}) \end{align}$

　各種エントロピーは非負の値なので、式(3)(4)から、結合エントロピーは各事象のエントロピーと条件付きエントロピー以上の値なのが分かります。

$\displaystyle \begin{aligned} H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{x}) \\ H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{y}) \\ H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{x} | \mathbf{y}) \\ H(\mathbf{x}, \mathbf{y}) &\geq H(\mathbf{y} | \mathbf{x}) \end{aligned}$

　以上で、各エントロピーとの関係が得られました。

可視化

　次は、2変数(2次元)の結合確率と自己情報量のヒートマップから結合エントロピーをグラフで確認します。この例では、2変数間に依存関係がないとします。依存関係を持つ場合の可視化については「条件付きエントロピーの定義 - からっぽのしょこ」を参照してください。

　利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)

　この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2を読み込む必要があります。
　また、ネイティブパイプ演算子|>を使っています。magrittrパッケージのパイプ演算子%>%に置き換えても処理できますが、その場合はmagrittrも読み込む必要があります。

インデックスの確認

　色々ややこしいので、インデックスとプロット位置の関係を確認しておきます。可視化自体には不要な処理です。

・作図コード(クリックで展開)

　事象系 $\mathbf{x}, \mathbf{y}$ それぞれの事象数を指定して、インデックスを作成します。

# 事象のインデックス(事象数)を指定
x_idx <- 1:4
y_idx <- 1:6

# xのインデックスを作成
index_x_df <- tibble::tibble(
  i   = x_idx, # xのインデックス
  p_x = 1 / length(x_idx), # xの確率:(等確率)
  I_x = - log2(p_x),       # xの自己情報量
  var_label = paste0("x[", i, "]") # 変数の式
)

# yのインデックスを作成
index_y_df <- tibble::tibble(
  j   = y_idx, # yのインデックス
  p_y = 1 / length(y_idx), # yの確率:(等確率)
  I_y = - log2(p_y),       # yの自己情報量
  var_label = paste0("y[", j, "]") # 変数の式
)
index_x_df; index_y_df

## # A tibble: 4 × 4
##       i   p_x   I_x var_label
##   <int> <dbl> <dbl> <chr>    
## 1     1  0.25     2 x[1]     
## 2     2  0.25     2 x[2]     
## 3     3  0.25     2 x[3]     
## 4     4  0.25     2 x[4]
## # A tibble: 6 × 4
##       j   p_y   I_y var_label
##   <int> <dbl> <dbl> <chr>    
## 1     1 0.167  2.58 y[1]     
## 2     2 0.167  2.58 y[2]     
## 3     3 0.167  2.58 y[3]     
## 4     4 0.167  2.58 y[4]     
## 5     5 0.167  2.58 y[5]     
## 6     6 0.167  2.58 y[6]

　1:事象数の形式で事象の数を指定して、1から事象数までの整数を作成しインデックス*_idxとします。
　インデックスをデータフレームに格納して、変数ラベルを作成します。数式を描画する場合はexpression()の記法を使います。下付き文字は変数[添字]の形式です。
　確認のため、事象ごとに一様な確率を作成して、自己情報量を計算していますが、ここでは使いません。

　インデックス $i$ について横方向に並べたグラフを作成します。

# xのインデックスを作図
index_x_graph <- ggplot() + 
  geom_tile(data = index_x_df, 
            mapping = aes(x = i, y = 0), 
            fill = "white", color = "black") + # インデックスのマップ
  geom_text(data = index_x_df, 
            mapping = aes(x = i, y = 0, label = var_label), parse = TRUE) + # 変数の式
  scale_x_continuous(breaks = x_idx, limits = c(0.5, max(x_idx)+0.5)) + # xのインデックス
  scale_y_continuous(breaks = 0, limits = c(-0.6, 0.6), labels = NULL) + # yのインデックス
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "index", 
       subtitle = expression(x == (list(x[1], cdots, x[n]))), 
       x = "i", y = "")
index_x_graph

　geom_tile()でヒートマップ用の事象セル、geom_text()で事象ラベルを描画します。数式を描画する場合はparse引数にTRUEを指定します。
　scale_*_continuous()のlimits引数に描画範囲の最小値と最大値を指定することでプロット領域のサイズを調整できます。

　同様に、インデックス $j$ について縦方向に並べたグラフを作成します。

# yのインデックスを作図
index_y_graph <- ggplot() + 
  geom_tile(data = index_y_df, 
            mapping = aes(x = 0, y = j), 
            fill = "white", color = "black") + # インデックスのマップ
  geom_text(data = index_y_df, 
            mapping = aes(x = 0, y = j, label = var_label), parse = TRUE) + # 変数の式
  scale_x_continuous(breaks = 0, limits = c(-0.6, 0.6), labels = NULL) + # xのインデックス
  scale_y_continuous(breaks = y_idx, limits = c(0.5, max(y_idx)+0.5)) + # yのインデックス
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "index", 
       subtitle = expression(y == (list(y[1], cdots, y[m]))), 
       x = "", y = "j")
index_y_graph

　2種類のインデックスの組み合わせを作成します。

# x・yのインデックスを作成
index_xy_df <- tidyr::expand_grid(
  i = x_idx, # xのインデックス
  j = y_idx  # yのインデックス
) |> # 全ての組み合わせを作成
  dplyr::mutate(
    p_x = 1 / length(unique(i)), # xの確率:(等確率)
    p_y = 1 / length(unique(j)), # yの確率:(等確率)
    p_xy = p_x * p_y,    # xyの結合確率
    I_xy = - log2(p_xy), # xyの自己情報量
    var_label = paste0("list(x[", i, "], y[", j, "])") # 変数の式
  )
index_xy_df

## # A tibble: 24 × 7
##        i     j   p_x   p_y   p_xy  I_xy var_label       
##    <int> <int> <dbl> <dbl>  <dbl> <dbl> <chr>           
##  1     1     1  0.25 0.167 0.0417  4.58 list(x[1], y[1])
##  2     1     2  0.25 0.167 0.0417  4.58 list(x[1], y[2])
##  3     1     3  0.25 0.167 0.0417  4.58 list(x[1], y[3])
##  4     1     4  0.25 0.167 0.0417  4.58 list(x[1], y[4])
##  5     1     5  0.25 0.167 0.0417  4.58 list(x[1], y[5])
##  6     1     6  0.25 0.167 0.0417  4.58 list(x[1], y[6])
##  7     2     1  0.25 0.167 0.0417  4.58 list(x[2], y[1])
##  8     2     2  0.25 0.167 0.0417  4.58 list(x[2], y[2])
##  9     2     3  0.25 0.167 0.0417  4.58 list(x[2], y[3])
## 10     2     4  0.25 0.167 0.0417  4.58 list(x[2], y[4])
## # … with 14 more rows

　インデックス $i, j$ を格納して、expand_grid()で全ての組み合わせを作成します。
　組み合わせごとに変数ラベルを作成します。複数の変数はlist(変数1, 変数2)の形式です。

　2種類のインデックス $i, j$ を横・縦方向に並べたグラフを作成します。

# x・yのインデックスの組み合わせを作図
index_xy_graph <- ggplot() + 
  geom_tile(data = index_xy_df, 
            mapping = aes(x = i, y = j), 
            fill = "white", color = "black") + # インデックスのマップ
  geom_text(data = index_xy_df, 
            mapping = aes(x = i, y = j, label = var_label), parse = TRUE) + # 変数の式
  scale_x_continuous(breaks = x_idx, limits = c(0.5, max(x_idx)+0.5)) + # xのインデックス
  scale_y_continuous(breaks = y_idx, limits = c(0.5, max(y_idx)+0.5)) + # yのインデックス
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "joint index", 
       subtitle = expression(list(x == (list(x[1], cdots, x[n])), y == (list(y[1], cdots, y[m])))), 
       x = "i", y = "j")
index_xy_graph

　先ほどと同様に作図します。

　3つのグラフを並べて描画します。

# グラフサイズ用の値を設定
x_size <- length(x_idx) - 1
y_size <- length(y_idx) - 1

# 並べて描画
index_joint_graph <- patchwork::wrap_plots(
  patchwork::plot_spacer(), index_x_graph, 
  index_y_graph, index_xy_graph, 
  nrow = 2, ncol = 2, 
  widths = c(1, x_size), heights = c(1, y_size)
)
index_joint_graph

　patchworkパッケージのwrap_plots()を使って、複数のグラフを並べて描画します。グラフを配置しない場所にはplot_spacer()を指定します。
　(この例では、アスペクト比を1にしているので、縦・横のサイズは事象数に依存しますが、(たぶんタイトルなどの領域の関係で)事象数のままだと微妙にズレたので、*_sizeの作成の際に値を調整しています。画像サイズの影響も受けます。)

　 $\mathbf{x}$ のインデックス $i$ について横方向、 $\mathbf{y}$ のインデックス $j$ について縦方向に並べて配置すると、格子点が全ての組み合わせになるのが分かります。
　セルの縦または横サイズを確率の値に、セルの色を確率またはエントロピーの値に対応させて可視化します。

結合確率のヒートマップ

　では、結合確率をヒートマップで可視化します。

・作図コード(クリックで展開)

　事象数を指定して、2つの事象系のインデックスと確率を作成します。

# 事象のインデックス(事象数)を指定
x_idx <- 1:4
y_idx <- 1:6

# 事象の確率を一様に作成
#p_x_vec <- rep(1 / length(x_idx), times = length(x_idx))
#p_y_vec <- rep(1 / length(y_idx), times = length(y_idx))

# 事象の確率をランダムに作成
p_x_vec <- MCMCpack::rdirichlet(n = 1, alpha = rep(1, times = length(x_idx))) |> 
  as.vector()
p_y_vec <- MCMCpack::rdirichlet(n = 1, alpha = rep(1, times = length(y_idx))) |> 
  as.vector()
p_x_vec; p_y_vec

## [1] 0.24374493 0.14587080 0.57909386 0.03129041
## [1] 0.20841422 0.34887220 0.01521039 0.18937211 0.06345578 0.17467530

　1:事象数の形で事象の数を指定して、インデックスを作成します。

　各事象の確率を等しく設定する場合は、(事象の数をnとして)n分の1の値をn個作成します。
　各事象の確率をランダムに設定する場合は、ディリクレ分布に従う乱数(総和が1になる非負の値)を使います。ディリクレ乱数は、MCMCpackパッケージのrdirichlet()で生成できます。サンプルサイズの引数nに1、パラメータの引数alphaにn個の値を指定します。

　事象 $x_i, y_j$ の値は可視化に影響しないので、ここでは省略します。

　確率に応じてプロットするための座標を作成します。

# 座標を作成
coord_x_vec <- c(0, cumsum(p_x_vec)[-length(p_x_vec)]) + 0.5*p_x_vec
coord_y_vec <- c(0, cumsum(p_y_vec)[-length(p_y_vec)]) + 0.5*p_y_vec
coord_x_vec; coord_y_vec

## [1] 0.1218725 0.3166803 0.6791627 0.9843548
## [1] 0.1042071 0.3828503 0.5648916 0.6671829 0.7935968 0.9126623

　各セルの縦または横サイズを確率に合わせます。そのため、各セルの中心の座標を「累積確率」+「その事象の確率の半分」の値とします。

　各事象の自己情報量を計算します。

# xの自己情報量を計算
entropy_x_df <- tibble::tibble(
  i   = x_idx, # xのインデックス
  p_x = p_x_vec,     # xの確率
  I_x = - log2(p_x), # xの自己情報量
  coord_x = coord_x_vec, # x軸の値
  p_label = paste0("p(x[", i, "])") # 確率の式
)
# yの自己情報量を計算
entropy_y_df <- tibble::tibble(
  j   = y_idx, # yのインデックス
  p_y = p_y_vec,     # yの確率
  I_y = - log2(p_y), # yの自己情報量
  coord_y = coord_y_vec, # y軸の値
  p_label = paste0("p(y[", j, "])") # 確率の式
)
entropy_x_df; entropy_y_df

## # A tibble: 4 × 5
##       i    p_x   I_x coord_x p_label
##   <int>  <dbl> <dbl>   <dbl> <chr>  
## 1     1 0.244  2.04    0.122 p(x[1])
## 2     2 0.146  2.78    0.317 p(x[2])
## 3     3 0.579  0.788   0.679 p(x[3])
## 4     4 0.0313 5.00    0.984 p(x[4])
## # A tibble: 6 × 5
##       j    p_y   I_y coord_y p_label
##   <int>  <dbl> <dbl>   <dbl> <chr>  
## 1     1 0.208   2.26   0.104 p(y[1])
## 2     2 0.349   1.52   0.383 p(y[2])
## 3     3 0.0152  6.04   0.565 p(y[3])
## 4     4 0.189   2.40   0.667 p(y[4])
## 5     5 0.0635  3.98   0.794 p(y[5])
## 6     6 0.175   2.52   0.913 p(y[6])

　インデックス $i$ と確率 $p(x_i)$ をデータフレームに格納して、自己情報量 $I(x_i) = - \log_2 p(x_i)$ を計算します。 $y_j$ についても同様です。
　座標列とラベル列も作成します。

　独立した2事象の自己情報量を計算します。

# x・yの自己情報量を計算
entropy_xy_df <- tidyr::expand_grid(
  i = x_idx, # xのインデックス
  j = y_idx  # yのインデックス
) |> # 全ての組み合わせを作成
  dplyr::mutate(
    p_x = p_x_vec[i], # xの確率
    p_y = p_y_vec[j], # yの確率
    coord_x = coord_x_vec[i], # x軸の値
    coord_y = coord_y_vec[j], # y軸の値
    p_xy = p_x * p_y,    # xyの結合確率
    I_xy = - log2(p_xy), # xyの自己情報量
    p_label = paste0("p(x[", i, "], y[", j, "])") # 確率の式
  )
entropy_xy_df

## # A tibble: 24 × 9
##        i     j   p_x    p_y coord_x coord_y    p_xy  I_xy p_label      
##    <int> <int> <dbl>  <dbl>   <dbl>   <dbl>   <dbl> <dbl> <chr>        
##  1     1     1 0.244 0.208    0.122   0.104 0.0508   4.30 p(x[1], y[1])
##  2     1     2 0.244 0.349    0.122   0.383 0.0850   3.56 p(x[1], y[2])
##  3     1     3 0.244 0.0152   0.122   0.565 0.00371  8.08 p(x[1], y[3])
##  4     1     4 0.244 0.189    0.122   0.667 0.0462   4.44 p(x[1], y[4])
##  5     1     5 0.244 0.0635   0.122   0.794 0.0155   6.01 p(x[1], y[5])
##  6     1     6 0.244 0.175    0.122   0.913 0.0426   4.55 p(x[1], y[6])
##  7     2     1 0.146 0.208    0.317   0.104 0.0304   5.04 p(x[2], y[1])
##  8     2     2 0.146 0.349    0.317   0.383 0.0509   4.30 p(x[2], y[2])
##  9     2     3 0.146 0.0152   0.317   0.565 0.00222  8.82 p(x[2], y[3])
## 10     2     4 0.146 0.189    0.317   0.667 0.0276   5.18 p(x[2], y[4])
## # … with 14 more rows

　インデックス $i, j$ を格納して、expand_grid()で全ての組み合わせを作成します。
　 $i, j$ の値を用いて、対応する確率と座標を取り出して格納します。
　結合確率 $p(x_i, y_j) = p(x_i) p(y_j)$ と自己情報量 $I(x_i, y_j) = - \log_2 p(x_i, y_j)$ を計算します。

　確率分布 $p(\mathbf{x})$ のヒートマップを作成します。

# 確率(グラデーション)の最大値を設定
p_max <- c(entropy_x_df[["p_x"]], entropy_y_df[["p_y"]], entropy_xy_df[["p_xy"]]) |> 
  max() |> 
  (\(.){ceiling(. * 10) / 10})() # 小数点第二位で切り上げ

# xの確率を作図
prob_x_graph <- ggplot() + 
  geom_tile(data = entropy_x_df, 
            mapping = aes(x = coord_x, y = 0, width = p_x, height = 1, fill = p_x), 
            color = "black", alpha = 0.8) + # 確率
  geom_text(data = entropy_x_df, 
            mapping = aes(x = coord_x, y = 0, label = p_label), parse = TRUE) + # 確率の式
  scale_x_continuous(breaks = coord_x_vec, labels = x_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_x_vec)), digits = 2), 
                                         name = expression(sum(p(x[i]), i)))) + # xのインデックスと累積確率
  scale_y_continuous(breaks = 0, labels = NULL) + # yのインデックス
  scale_fill_viridis_c(limits = c(0, p_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "probability", 
       subtitle = expression(list(x == (list(x[1], cdots, x[n])), p(x))), 
       fill = "p", 
       x = "i", y = "")
prob_x_graph

　他のヒートマップとグラデーションの基準を統一するために、確率の最大値をp_maxとしておきます。

　geom_tile()でヒートマップを描画します。x引数に座標列coord_x、width引数に確率列p_xを指定して、各セルの横サイズを確率に対応付けます。y引数の値は何でもよく、height引数に1を指定して、ヒートマップ全体の面積が1×1になるようにします。また、fill引数に確率列を指定して、各セルの色を確率に対応付けます。
　scale_*_continuous()で軸目盛を設定します。breaks引数に確率に応じた座標coord_x_vec、labels引数にインデックスx_idxを指定して、主軸(第1軸)目盛が各セルの中心に対応するようにインデックスを表示します。sec.axis引数にsec_axis()を使って、第2軸が各セルの端に対応するように累積確率を表示します。

　同様に、確率分布 $p(\mathbf{y})$ のヒートマップを作成します。

# yの確率を作図
prob_y_graph <- ggplot() + 
  geom_tile(data = entropy_y_df, 
            mapping = aes(x = 0, y = coord_y, width = 1, height = p_y, fill = p_y), 
            color = "black", alpha = 0.8) + # 確率
  geom_text(data = entropy_y_df, 
            mapping = aes(x = 0, y = coord_y, label = p_label), parse = TRUE) + # 確率の式
  scale_x_continuous(breaks = 0, labels = NULL) + # xのインデックス
  scale_y_continuous(breaks = coord_y_vec, labels = y_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_y_vec)), digits = 2), 
                                         name = expression(sum(p(y[j]), j)))) + # yのインデックスと累積確率
  scale_fill_viridis_c(limits = c(0, p_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "probability", 
       subtitle = expression(list(y == (list(y[1], cdots, y[m])), p(y))), 
       fill = "p", 
       x = "", y = "j")
prob_y_graph

　 $p(\mathbf{x}, \mathbf{y})$ のヒートマップを作成します。

# x・yの結合確率を作図
prob_xy_graph <- ggplot() + 
  geom_tile(data = entropy_xy_df, 
            mapping = aes(x = coord_x, y = coord_y, width = p_x, height = p_y, fill = p_xy), 
            color = "black", alpha = 0.8) + # 確率
  geom_text(data = entropy_xy_df, 
            mapping = aes(x = coord_x, y = coord_y, label = p_label), parse = TRUE, 
            angle = 0) + # 確率の式
  scale_x_continuous(breaks = coord_x_vec, labels = x_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_x_vec)), digits = 2), 
                                         name = expression(sum(p(x[i]), i)))) + # xのインデックスと累積確率
  scale_y_continuous(breaks = coord_y_vec, labels = y_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_y_vec)), digits = 2), 
                                         name = expression(sum(p(y[j]), j)))) + # yのインデックスと累積確率
  scale_fill_viridis_c(limits = c(0, p_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "joint probability", 
       subtitle = expression(p(x[i], y[j]) == p(x[i]) * p(y[j])), 
       fill = "p", 
       x = "i", y = "j")
prob_xy_graph

　3つのグラフを並べて描画します。

# 並べて描画
prob_joint_graph <- patchwork::wrap_plots(
  patchwork::plot_spacer(), prob_x_graph, 
  prob_y_graph, prob_xy_graph, 
  nrow = 2, ncol = 2, 
  guides = "collect"
)
prob_joint_graph

　wrap_plots()でグラフを並べて描画します。グラフを表示しない位置にはplot_spacer()を指定します。

　アスペクト比を1にしてセルの縦または横サイズを確率に合わせることで、全体の面積に占めるセル面積の割合が各事象の確率に対応します。セルが大きいほど確率が大きく、セルが小さいほど確率も小さいことを表します。

　 $p(x_i)$ の値(横方向のサイズ)と $p(y_j)$ の値(縦方向のサイズ)により $p(x_i, y_j)$ の値(セルサイズ)が決まり、2事象に関する総和が1(全体の面積が1×1) $\sum_{i=1}^n \sum_{j=1}^m = p(x_i, y_j) = 1$ になるのが分かります。

自己情報量のヒートマップ

　続いて、自己情報量をヒートマップで可視化して、結合エントロピーを確認します。

・作図コード(クリックで展開)

　 $\mathbf{x}$ の各事象の自己情報量のヒートマップを作成します。

# 自己情報量(グラデーション)の最大値を設定
I_max <- c(entropy_x_df[["I_x"]], entropy_y_df[["I_y"]], entropy_xy_df[["I_xy"]]) |> 
  #(\(.){.[. < Inf]})() |> # 確率が0の事象を除く
  max() |> 
  ceiling() # 小数点以下を切り上げ

# 平均情報量を計算
H_x <- sum(entropy_x_df[["p_x"]] * entropy_x_df[["I_x"]], na.rm = TRUE)

# ラベル用の文字列を作成
fnc_label <- paste0(
  "list(", 
  "I(x[i]) == - log[2]*p(x[i])", ", ", 
  "H(x) == sum(p(x[i]) * I(x[i]), i==1, n)", 
  ")"
)
H_x_label <- paste0("H(x) == ", round(H_x, digits = 3))

# xの自己情報量を作図
entropy_x_graph <- ggplot() + 
  geom_tile(data = entropy_x_df, 
            mapping = aes(x = coord_x, y = 0, width = p_x, height = 1, fill = I_x), 
            color = "black", alpha = 0.8) + # 自己情報量
  geom_text(data = entropy_x_df, 
            mapping = aes(x = coord_x, y = 0, label = p_label), parse = TRUE) + # 確率の式
  geom_label(mapping = aes(x = -Inf, y = Inf), label = H_x_label, hjust = 0, vjust = 1, parse = TRUE, 
             alpha = 0.5) + # エントロピーの値
  scale_x_continuous(breaks = coord_x_vec, labels = x_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_x_vec)), digits = 2), 
                                         name = expression(sum(p(x[i]), i)))) + # xのインデックスと累積確率
  scale_y_continuous(breaks = 0, labels = NULL) + # yのインデックス
  scale_fill_viridis_c(limits = c(0, I_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "self-information", 
       subtitle = parse(text = fnc_label), 
       fill = "I", 
       x = "i", y = "")
entropy_x_graph

　自己情報量の最大値をI_maxとしておき、「結合確率のヒートマップ」のときと同様に作図します。
　塗りつぶし色(z軸の値)を自己情報量列I_xとします。

　エントロピー(平均情報量) $H(\mathbf{x}) = - \sum_{i=1}^n p(x_i) I(x_i)$ を計算して、geom_label()でラベルとして表示します。

　同様に、 $\mathbf{y}$ の各事象の自己情報量のヒートマップを作成します。

# 平均情報量を計算
H_y <- sum(entropy_y_df[["p_y"]] * entropy_y_df[["I_y"]], na.rm = TRUE)

# ラベル用の文字列を作成
fnc_label <- paste0(
  "list(", 
  "I(y[j]) == - log[2]*p(y[j])", ", ", 
  "H(y) == sum(p(y[j]) * I(y[j]), j==1, m)", 
  ")"
)
H_y_label <- paste0("H(y) == ", round(H_y, digits = 3))

# yの自己情報量を作図
entropy_y_graph <- ggplot() + 
  geom_tile(data = entropy_y_df, 
            mapping = aes(x = 0, y = coord_y, width = 1, height = p_y, fill = I_y), 
            color = "black", alpha = 0.8) + # 自己情報量
  geom_text(data = entropy_y_df, 
            mapping = aes(x = 0, y = coord_y, label = p_label), parse = TRUE) + # 確率の式
  geom_label(mapping = aes(x = -Inf, y = Inf), label = H_y_label, hjust = 0, vjust = 1, parse = TRUE, 
             alpha = 0.5) + # エントロピーの値
  scale_x_continuous(breaks = 0, labels = NULL) + # xのインデックス
  scale_y_continuous(breaks = coord_y_vec, labels = y_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_y_vec)), digits = 2), 
                                         name = expression(sum(p(y[j]), j)))) + # yのインデックスと累積確率
  scale_fill_viridis_c(limits = c(0, I_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "self-information", 
       subtitle = parse(text = fnc_label), 
       fill = "I", 
       x = "", y = "j")
entropy_y_graph

　 $\mathbf{x}, \mathbf{y}$ の事象の組み合わせごとの自己情報量のヒートマップを作成します。

# 平均情報量を計算
H_xy <- sum(entropy_xy_df[["p_xy"]] * entropy_xy_df[["I_xy"]], na.rm = TRUE)

# ラベル用の文字列を作成
fnc_label <- paste0(
  "list(", 
  "I(x[i], y[j]) == - log[2]*p(x[i], y[j])", ", ", 
  "H(x, y) == sum({}, i==1, n) * sum({}, j==1, m) * p(x[i], y[j]) * I(x[i], y[j])", 
  ")"
)
H_xy_label <- paste0("H(x, y) == ", round(H_xy, digits = 3))

# x・yの自己情報量を作図
entropy_xy_graph <- ggplot() + 
  geom_tile(data = entropy_xy_df, 
            mapping = aes(x = coord_x, y = coord_y, width = p_x, height = p_y, fill = I_xy), 
            color = "black", alpha = 0.8) + # 自己情報量
  geom_text(data = entropy_xy_df, 
            mapping = aes(x = coord_x, y = coord_y, label = p_label), parse = TRUE, 
            angle = 0) + # 確率の式
  geom_label(mapping = aes(x = -Inf, y = Inf), label = H_xy_label, hjust = 0, vjust = 1, parse = TRUE, 
             alpha = 0.5) + # エントロピーの値
  scale_x_continuous(breaks = coord_x_vec, labels = x_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_x_vec)), digits = 2), 
                                         name = expression(sum(p(x[i]), i)))) + # xのインデックスと累積確率
  scale_y_continuous(breaks = coord_y_vec, labels = y_idx, 
                     sec.axis = sec_axis(trans = ~., breaks = round(cumsum(c(0, p_y_vec)), digits = 2), 
                                         name = expression(sum(p(y[j]), j)))) + # yのインデックスと累積確率
  scale_fill_viridis_c(limits = c(0, I_max)) + # グラデーション
  coord_fixed(ratio = 1) + # アスペクト比
  theme(panel.grid.minor = element_blank()) + # 図の体裁
  labs(title = "joint information", 
       subtitle = parse(text = fnc_label), 
       fill = "I", 
       x = "i", y = "j")
entropy_xy_graph

　結合エントロピー(多変数の平均情報量) $H(\mathbf{x}, \mathbf{y}) = - \sum_{i=1}^n \sum_{j=1}^m p(x_i, y_j) I(x_i, y_j)$ を計算して、ラベルとして表示します。

　3つのグラフを並べて描画します。

# 並べて描画
entropy_joint_graph <- patchwork::wrap_plots(
  patchwork::plot_spacer(), entropy_x_graph, 
  entropy_y_graph, entropy_xy_graph, 
  nrow = 2, ncol = 2, 
  guides = "collect"
)
entropy_joint_graph

　確率(セルサイズ)が小さいほど、自己情報量が大きくなります。各事象の確率と自己情報量の関係については「自己情報量の定義 - からっぽのしょこ」を参照してください。
　事象(セル)ごとの自己情報量(色を高さとみなして)の期待値(面積の割合と高さの加重平均)がエントロピーや結合エントロピーの値です。結合エントロピーは、1つの事象の確率が1で他の事象の確率が全て0のとき最小値の $H_{\min}(\mathbf{x}, \mathbf{y}) = 0$ 、全ての事象の確率が等しい $p(x_i, y_j) = \frac{1}{n m}$ とき最大値の $H_{\max}(\mathbf{x}, \mathbf{y}) = - \log_2 (n m)$ になります。全ての事象の確率と自己情報量の関係については「平均情報量の最大値(最大エントロピー)の導出 - からっぽのしょこ」を参照してください。