Rで生成
# パッケージ読み込み
library(tidyverse)
# 乱数固定
set.seed(20230228)
乱数生成のためのサンプルサイズ設定とパラメータ設定.
# サンプルサイズ
N <- 10000
# ディリクレ分布のパラメータ設定
alpha <- c(1, 2, 3)
ガンマ分布から乱数を生成
仕組みを分かりやすくするために, 無駄にループを回す, Rとしては非効率な書き方になっているので注意(後述)
k <- length(alpha)
z <- matrix(numeric(k*N), ncol=k, byrow=TRUE)
q_gamma <- matrix(numeric(k*N), ncol=k, byrow=TRUE)
# Step 1 : ガンマ乱数生成
for(i in 1:k){
z[,i] <- rgamma(n=N, alpha[i])
}
# Step 2 : 正規化
for(j in 1:N){
q_gamma[j,] <- z[j,]/sum(z[j,])
}
{MCMCpack}パッケージの関数を使用して生成
q_mcmcpack <- MCMCpack::rdirichlet(n = N, alpha = alpha)
生成したディリクレ乱数を可視化したいが, 自分では作画できないので, 以下の記事のコードをほとんどそのまま使用させていただきました. 感謝. ディリクレ分布の乱数生成
www.anarchive-beta.com
自分でガンマ乱数から生成したもののヒートマップを描く.
# 軸目盛の位置を指定
axis_vals <- seq(from = 0, to = 1, by = 0.1)
# 枠線用の値を作成
ternary_axis_df <- tibble::tibble(
y_1_start = c(0.5, 0, 1), # 始点のx軸の値
y_2_start = c(0.5*sqrt(3), 0, 0), # 始点のy軸の値
y_1_end = c(0, 1, 0.5), # 終点のx軸の値
y_2_end = c(0, 0, 0.5*sqrt(3)), # 終点のy軸の値
axis = c("x_1", "x_2", "x_3") # 元の軸
)
# グリッド線用の値を作成
ternary_grid_df <- tibble::tibble(
y_1_start = c(
0.5 * axis_vals,
axis_vals,
0.5 * axis_vals + 0.5
), # 始点のx軸の値
y_2_start = c(
sqrt(3) * 0.5 * axis_vals,
rep(0, times = length(axis_vals)),
sqrt(3) * 0.5 * (1 - axis_vals)
), # 始点のy軸の値
y_1_end = c(
axis_vals,
0.5 * axis_vals + 0.5,
0.5 * rev(axis_vals)
), # 終点のx軸の値
y_2_end = c(
rep(0, times = length(axis_vals)),
sqrt(3) * 0.5 * (1 - axis_vals),
sqrt(3) * 0.5 * rev(axis_vals)
), # 終点のy軸の値
axis = c("x_1", "x_2", "x_3") |>
rep(each = length(axis_vals)) # 元の軸
)
# 軸ラベル用の値を作成
ternary_axislabel_df <- tibble::tibble(
y_1 = c(0.25, 0.5, 0.75), # x軸の値
y_2 = c(0.25*sqrt(3), 0, 0.25*sqrt(3)), # y軸の値
label = c("alpha[1]", "alpha[2]", "alpha[3]"), # 軸ラベル
h = c(3, 0.5, -2), # 水平方向の調整用の値
v = c(0.5, 3, 0.5), # 垂直方向の調整用の値
axis = c("x_1", "x_2", "x_3") # 元の軸
)
# 軸目盛ラベル用の値を作成
ternary_ticklabel_df <- tibble::tibble(
y_1 = c(
0.5 * axis_vals,
axis_vals,
0.5 * axis_vals + 0.5
), # x軸の値
y_2 = c(
sqrt(3) * 0.5 * axis_vals,
rep(0, times = length(axis_vals)),
sqrt(3) * 0.5 * (1 - axis_vals)
), # y軸の値
label = c(
rev(axis_vals),
axis_vals,
rev(axis_vals)
), # 軸目盛ラベル
h = c(
rep(1.5, times = length(axis_vals)),
rep(1.5, times = length(axis_vals)),
rep(-0.5, times = length(axis_vals))
), # 水平方向の調整用の値
v = c(
rep(0.5, times = length(axis_vals)),
rep(0.5, times = length(axis_vals)),
rep(0.5, times = length(axis_vals))
), # 垂直方向の調整用の値
angle = c(
rep(-60, times = length(axis_vals)),
rep(60, times = length(axis_vals)),
rep(0, times = length(axis_vals))
), # ラベルの表示角度
axis = c("x_1", "x_2", "x_3") |>
rep(each = length(axis_vals)) # 元の軸
)
# サンプルを三角座標に変換して格納(q_gamma)
data_df_q_gamma <- tibble::tibble(
y_1 = q_gamma[, 2] + 0.5 * q_gamma[, 3], # 三角座標のx軸の値
y_2 = sqrt(3) * 0.5 * q_gamma[, 3] # 三角座標のy軸の値
)
# 三角座標の値を作成
y_1_vals <- seq(from = 0, to = 1, length.out = 301)
y_2_vals <- seq(from = 0, to = 0.5*sqrt(3), length.out = 300)
# 格子点を作成
y_mat <- tidyr::expand_grid(
y_1 = y_1_vals,
y_2 = y_2_vals
) |> # 格子点を作成
as.matrix() # マトリクスに変換
# 3次元変数に変換
phi_mat <- tibble::tibble(
phi_2 = y_mat[, 1] - y_mat[, 2] / sqrt(3),
phi_3 = 2 * y_mat[, 2] / sqrt(3)
) |> # 元の座標に変換
dplyr::mutate(
phi_2 = dplyr::if_else(phi_2 >= 0 & phi_2 <= 1, true = phi_2, false = as.numeric(NA)),
phi_3 = dplyr::if_else(phi_3 >= 0 & phi_3 <= 1 & !is.na(phi_2), true = phi_3, false = as.numeric(NA)),
phi_1 = 1 - phi_2 - phi_3,
phi_1 = dplyr::if_else(phi_1 >= 0 & phi_1 <= 1, true = phi_1, false = as.numeric(NA))
) |> # 範囲外の値をNAに置換
dplyr::select(phi_1, phi_2, phi_3) |> # 順番を変更
as.matrix() # マトリクスに変換
# パラメータラベル用の文字列を作成
param_text <- paste0(
"list(",
"alpha==(list(", paste0(alpha, collapse = ", "), "))",
", N==", N,
")"
)
# サンプルの度数のヒートマップを作成
ggplot() +
geom_segment(data = ternary_axis_df,
mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end),
color = "gray50") + # 三角図の枠線
geom_segment(data = ternary_grid_df,
mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end),
color = "gray50", linetype = "dashed") + # 三角図のグリッド線
geom_text(data = ternary_ticklabel_df,
mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v, angle = angle)) + # 三角図の軸目盛ラベル
geom_text(data = ternary_axislabel_df,
mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v),
parse = TRUE, size = 6) + # 三角図の軸ラベル
geom_bin_2d(data = data_df_q_gamma,
mapping = aes(x = y_1, y = y_2, fill = ..count..),
alpha = 0.8) + # サンプル
scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
scale_x_continuous(breaks = c(0, 0.5, 1), labels = NULL) + # x軸
scale_y_continuous(breaks = c(0, 0.25*sqrt(3), 0.5*sqrt(3)), labels = NULL) + # y軸
coord_fixed(ratio = 1, clip = "off") + # アスペクト比
theme(axis.ticks = element_blank(),
panel.grid.minor = element_blank()) + # 図の体裁
labs(title = "Dirichlet Distribution (Gamma RVs method)",
subtitle = parse(text = param_text),
fill = "frequency",
x = "", y = "") +
theme_bw()
# サンプルを三角座標に変換して格納(q_mcmcpack)
data_df_q_mcmcpack <- tibble::tibble(
y_1 = q_mcmcpack[, 2] + 0.5 * q_mcmcpack[, 3], # 三角座標のx軸の値
y_2 = sqrt(3) * 0.5 * q_mcmcpack[, 3] # 三角座標のy軸の値
)
# サンプルの度数のヒートマップを作成
ggplot() +
geom_segment(data = ternary_axis_df,
mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end),
color = "gray50") + # 三角図の枠線
geom_segment(data = ternary_grid_df,
mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end),
color = "gray50", linetype = "dashed") + # 三角図のグリッド線
geom_text(data = ternary_ticklabel_df,
mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v, angle = angle)) + # 三角図の軸目盛ラベル
geom_text(data = ternary_axislabel_df,
mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v),
parse = TRUE, size = 6) + # 三角図の軸ラベル
geom_bin_2d(data = data_df_q_mcmcpack,
mapping = aes(x = y_1, y = y_2, fill = ..count..),
alpha = 0.8) + # サンプル
scale_fill_distiller(palette = "Spectral") + # 塗りつぶしの色
scale_x_continuous(breaks = c(0, 0.5, 1), labels = NULL) + # x軸
scale_y_continuous(breaks = c(0, 0.25*sqrt(3), 0.5*sqrt(3)), labels = NULL) + # y軸
coord_fixed(ratio = 1, clip = "off") + # アスペクト比
theme(axis.ticks = element_blank(),
panel.grid.minor = element_blank()) + # 図の体裁
labs(title = "Dirichlet Distribution (MCMCpack)",
subtitle = parse(text = param_text),
fill = "frequency",
x = "", y = "") +
theme_bw()
設定したパラメータから計算される平均
alpha/sum(alpha)
## [1] 0.1666667 0.3333333 0.5000000
自分でガンマ乱数から生成したものの平均値
colSums(q_gamma)/N
## [1] 0.1645176 0.3340670 0.5014154
{MCMCpack}パッケージの関数を使用して生成したものの平均値
colSums(q_mcmcpack)/N
## [1] 0.1665876 0.3357620 0.4976504
グラフ, 平均値の比較をみてみたところ, 合っていそう.
