第7章 パターンをとらえる指標化

図の右上のshowボタンを押すとRのコードが表示されます。

7.1 時間的なパターンをとらえる

7.1.1 じゃんけんマシンの分析

7.1.1.1

library(conflicted)
library(tidyverse)
library(patchwork)

# シード値を固定
set.seed(42)

# データ生成関数
generate_data <- function(transition_prob, n=1000) {
  initial_state <- sample.int(3, 1) # 初期状態をランダムに選択
  # 遷移確率に従って次の状態を選択
  states <- accumulate(
    .init = initial_state,
    .x = rep(1, n-1),
    .f = function(state, .) {
      sample.int(3, size = 1, prob = transition_prob[[state]])
      }
    )
  return(states)
}

# じゃんけんマシンAの遷移確率
transition_prob_A <- list(c(0.2, 0.5, 0.3), c(0.4, 0.1, 0.5), c(0.3, 0.6, 0.1))

# じゃんけんマシンBの遷移確率
transition_prob_B <- list(c(0.3, 0.4, 0.3), c(0.2, 0.3, 0.5), c(0.4, 0.2, 0.4))

# データ生成
data_A <- generate_data(transition_prob_A)
data_B <- generate_data(transition_prob_B)

# 最初の50回のじゃんけんの時系列を折れ線グラフで描画
df_long <- data.frame(Round=1:50, Machine_A=data_A[1:50], Machine_B=data_B[1:50]) |>
  pivot_longer(!Round)

p1 <- df_long |>
  ggplot(aes(x=Round, y=value, color=name)) +
  geom_line(alpha = 0.5) +
  geom_point(alpha = 0.5) +  
  scale_y_continuous(breaks=1:3, labels=c("グー", "チョキ", "パー")) +
  scale_color_discrete(labels = c("じゃんけんマシンA", "じゃんけんマシンB")) +
  labs(title="分析対象のジャンケン時系列（最初の50ラウンド）") +
  theme(
    axis.title = element_blank(),
    legend.title = element_blank(),
    legend.position = "bottom"
  )

# 各手の連続回数の平均値を計算する関数
calculate_run_mean <- function(data) {
  rle_data <- rle(data)  # 要素が連続する部分の値と長さを計算
  counts <- split(rle_data$lengths, rle_data$values)  # 各手の連続回数を取得
  return(sapply(counts, mean, na.rm=TRUE))  # 各手の連続回数の平均値を計算
}

# じゃんけんマシンAとBの各手の連続回数の平均値を計算
mean_run_A <- calculate_run_mean(data_A)
mean_run_B <- calculate_run_mean(data_B)

# 集団棒グラフで各手の出現割合を描画
p2 <- rbind(
  table(data_A) / length(data_A),
  table(data_B) / length(data_B)
  ) |>
  as.data.frame() |>
  mutate(machine = c("じゃんけんマシンA", "じゃんけんマシンB")) |>
  pivot_longer(!machine) |>
  ggplot(aes(x = name, y = value, fill = machine)) +
  geom_col(position = "dodge") +
  scale_x_discrete(labels = c("1" = "グー", "2" = "チョキ", "3" = "パー")) +
  theme(legend.position = "none", axis.title = element_blank()) +
  labs(title = "出現頻度")

# 集団棒グラフで各手の連続回数の平均値を描画
p3 <-rbind(mean_run_A, mean_run_B) |>
  as.data.frame() |>
  mutate(machine = c("じゃんけんマシンA", "じゃんけんマシンB")) |>
  pivot_longer(!machine) |>
  ggplot(aes(x = name, y = value, fill = machine)) +
  geom_col(position = "dodge") +
  scale_x_discrete(labels = c("1" = "グー", "2" = "チョキ", "3" = "パー")) +
  theme(legend.title = element_blank(), axis.title = element_blank()) +
  labs(title = "持続回数の平均値")

p1 / {p2 | p3}

7.1.1.2

ggraphで書き直したい

library(conflicted)
library(tidyverse)
library(igraph)
library(patchwork)
library(viridisLite)

# シード値を固定
set.seed(42)

# データ生成関数
generate_data <- function(transition_prob, n=1000) {
  initial_state <- sample.int(3, 1) # 初期状態をランダムに選択
  # 遷移確率に従って次の状態を選択
  states <- accumulate(
    .init = initial_state,
    .x = rep(1, n-1),
    .f = function(state, .) {
      sample.int(3, size = 1, prob = transition_prob[[state]])
      }
    )
  return(states)
}

# じゃんけんマシンAの遷移確率
transition_prob_A <- list(c(0.2, 0.5, 0.3), c(0.4, 0.1, 0.5), c(0.3, 0.6, 0.1))

# じゃんけんマシンBの遷移確率
transition_prob_B <- list(c(0.3, 0.4, 0.3), c(0.2, 0.3, 0.5), c(0.4, 0.2, 0.4))

# データ生成
data_A <- generate_data(transition_prob_A)
data_B <- generate_data(transition_prob_B)

# データから遷移確率を計算する関数
calculate_empirical_transition_prob <- function(data) {
  # 現在の状態と次の状態のペアを作成、遷移回数をカウント、遷移確率を計算
  transition_prob <- paste0(data[-length(data)], "_", data[-1]) |>
    table() |>
    prop.table() |>
    matrix(nrow = 3, byrow = TRUE)
  return(transition_prob)
}

# 実際のデータから遷移確率を計算
transition_prob_A <- calculate_empirical_transition_prob(data_A)
transition_prob_B <- calculate_empirical_transition_prob(data_B)

# ネットワーク生成関数
plot_transition_graph <- function(prob_matrix, title) {
  normalize_vector <- function(x) {  # ベクトルを標準化（最小値1、最大値10）
    min_x <- min(x)
    max_x <- max(x)
    normalized_x <- 1 + 99 * ((x - min_x) / (max_x - min_x))
    return(normalized_x)
    }

  my_cols <- turbo(100)
  
  combinations <- expand.grid(i = 1:3, j = 1:3) # 全ての組み合わせを生成
  g <- graph.empty(n=3, directed=TRUE) |>       # 空の有向グラフを作成
    add_edges(as.vector(t(combinations)))       # エッジを追加
  V(g)$name <- c("グー", "チョキ", "パー")      # ノードの名前を設定
  V(g)$size <- 40
  # エッジのラベルと太さを計算
  edge_labels <- sapply(
    1:nrow(combinations),
    function(x) sprintf("%.2f", prob_matrix[combinations[x, 1], combinations[x, 2]])
    )
  edge_widths <- sapply(
    1:nrow(combinations),
    function(x) prob_matrix[combinations[x, 1], combinations[x, 2]] * 20
    )
  edge_colors <- my_cols[normalize_vector(edge_widths)]
  E(g)$curved <- 0.5
  E(g)$label <- edge_labels
  E(g)$width <- edge_widths
  E(g)$color <- edge_colors
  g |>
    plot(main=title)  # ネットワークを描画
}

par(mfrow = c(1, 2), mar = c(0,1,1,1))
plot_transition_graph(transition_prob_A, "遷移確率（じゃんけんマシンA）")
plot_transition_graph(transition_prob_B, "遷移確率（じゃんけんマシンB）")

par(mfrow = c(1, 1), mar = c(5.1, 4.1, 4.1, 2.1))

7.1.2 自己相関で周期的なパターンを検出する

自己相関プロットの再現がイマイチ。Bartlett’s formulaを用いた95%信頼区間の表示が同じようにできない。

library(conflicted)
library(tidyverse)
library(patchwork)
library(forecast)

set.seed(42)
n <- 100
x <- seq(0, 100, length.out = n)

# データフレームにデータを格納
df <- data.frame(
  x = 1:n,
  y = sin(pi / 5 * x) + 0.8 * rnorm(n)
  )

# 元の時系列データをプロット
p1 <- df |>
  ggplot(aes(x = x, y = y)) +
  geom_line() +
  labs(
    x = expression(paste("時刻", italic(t))), 
    y = expression(italic(x(t))),
    title = "時系列データ"
    ) +
  theme(aspect.ratio = 1/3)

# 自己相関プロットを作成
p2 <- ggAcf(df$y, lag.max = 40) +
  labs(title = "自己相関プロット", x = "ラグ（ずれ幅）", y = "自己相関係数") +
  theme(aspect.ratio = 1/3)

# ラグ10のラグプロット
y <- df$y
y_lag10 <- dplyr::lead(y, 10)  # ラグ10を取る

# 相関係数を計算
corr <- cor(y, y_lag10, use = "pairwise.complete.obs")

# ラグプロットの作成
p3 <- data.frame(y, y_lag10) |>
  drop_na() |>
  ggplot(aes(x = y, y = y_lag10)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = "lm", color = "darkblue", fill = "darkblue") +
  annotate("text", x=-2.5, y=2, label=paste0("r = ", round(corr, 2))) +
  labs(x = expression(italic(x(t))), y = expression(italic(x(t+10))), title = "ラグ10での相関") +
  theme(aspect.ratio = 1)

(p1 + plot_spacer() + plot_layout(ncol = 2, widths = c(.7, .3))) / 
  (p2 + p3 + plot_layout(ncol = 2, widths = c(.7, .3)))

7.1.3 フーリエ変換で周波数の情報を取り出す

音声データの取り扱いが分からないので割愛

7.1.4 周波数の時間変化を見る

音声データの取り扱いが分からないので割愛

7.1.5 心電図データの分析

心電図データの取り扱いが分からないので割愛

7.2 空間データのパターンをとらえる

7.2.1 ボロノイ図による空間の分割

最大距離を制限した場合の各領域の面積の出し方が分からない

library(conflicted)
library(tidyverse)
library(deldir)
library(ggforce)
library(viridis)
library(patchwork)

set.seed(42)
df <- as.data.frame(
  rbind(
    matrix(c(0, 0, 0, 1, 1, 0, 1, 1), ncol = 2), #四隅
    matrix(runif(40), ncol = 2)
    )
  ) |>
  distinct()

# 最短距離
df2 <- df |>
  dist() |>
  as.matrix() |>
  as.data.frame() |>
  rowid_to_column() |>
  pivot_longer(!rowid) |>
  dplyr::filter(value > 0) |>
  slice_min(value, by = rowid)

df3 <- df2 |>
  left_join(
    df |>
      rowid_to_column(),
    by = join_by(rowid == rowid)
    ) |>
  left_join(
    df |>
      rowid_to_column() |>
      mutate(rowid = as.character(rowid)),
    by = join_by(name == rowid)
    ) |>
  mutate(
    x = V1.x, y = V2.x, xend = V1.y, yend = V2.y, .keep = "unused"
  )

res <- deldir(x = df3$x, y = df3$y)

df_area <- df3 |>
  mutate(area = res$summary$dir.area)

p1 <- df_area |>
  ggplot(aes(x = x, y = y, fill = area)) +
  geom_voronoi_tile(colour = "black") +
  geom_point() +
  geom_segment(
    aes(x = x, y = y, xend = xend, yend = yend),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "inches")
      ),
    color = "red"
    ) +
  scale_fill_viridis(option = "turbo") + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  theme(aspect.ratio = 1) +
  labs(title = "ボロノイ図", fill = "面積")

p2 <- df_area |>
  ggplot(aes(x = x, y = y, fill = area)) +
  geom_voronoi_tile(colour = "black", max.radius = 0.2) +
  geom_point() +
  scale_fill_viridis(option = "turbo") + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1))+
  theme(aspect.ratio = 1) +
  labs(title = "ボロノイ図（最大距離制限付き）", fill = "面積")

p3 <- df2 |>
  ggplot(aes(x = value)) +
  geom_histogram(bins = 10, fill = "blue", alpha = 0.7, color = "black") +
  xlab("最近傍点までの最短") +
  ylab("頻度")

p4 <- df_area |>
  ggplot(aes(x = area)) +
  geom_histogram(bins = 10, alpha = 0.7, color = "black") +
  xlab("各領域の面積") +
  ylab("頻度")

{p1 | p2 } / {p3 | p4}

7.2.2 カーネル密度推定による空間の利用密度の定量化

library(conflicted)
library(tidyverse)
library(deldir)
library(ggforce)
library(viridis)
library(patchwork)

set.seed(0)
df <- as.data.frame(
  rbind(
    matrix(c(0, 0, 0, 1, 1, 0, 1, 1), ncol = 2), #四隅
    matrix(runif(40), ncol = 2)
    )
  ) |>
  distinct()

# 新しいボロノイ図を計算
res <- deldir(x = df$V1, y = df$V2)

p1 <- df |>
  mutate(area = res$summary$dir.area) |>
  ggplot(aes(x = V1, y = V2, fill = area)) +
  geom_voronoi_tile(colour = "black") +
  geom_point() +
  scale_fill_viridis(option = "turbo") + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  theme(aspect.ratio = 1, legend.position = "none", axis.title = element_blank()) +
  labs(title = "ボロノイ図")

p2 <- df |>
  ggplot(aes(x = V1, y = V2)) +
  geom_density_2d_filled(h = 0.25) +
  geom_point(alpha = 0.5) +
  scale_fill_manual(values = turbo(7)) + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  theme(aspect.ratio = 1, legend.position = "none", axis.title = element_blank()) +
  labs(title = "KDEバンド幅=0.25")

p3 <- df |>
  ggplot(aes(x = V1, y = V2)) +
  geom_density_2d_filled(h = 0.5) +
  geom_point(alpha = 0.5) + 
  scale_fill_manual(values = turbo(8)) + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  theme(aspect.ratio = 1, legend.position = "none", axis.title = element_blank()) +
  labs(title = "KDEバンド幅=0.50")

p4 <- df |>
  ggplot(aes(x = V1, y = V2)) +
  geom_density_2d_filled(h = 0.75) +
  geom_point(alpha = 0.5) +
  scale_fill_manual(values = turbo(9)) + 
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  theme(aspect.ratio = 1, legend.position = "none", axis.title = element_blank()) +
  labs(title = "KDEバンド幅=0.75")


{p1|p2}/{p3|p4}

7.2.3 グレイレベル共起行列を用いた分析

画像の取り扱いがわからないので割愛

7.3 ネットワークのパターンをとらえる

7.3.1 様々な中心性指標

library(conflicted)
library(tidyverse)
library(igraph)
library(viridis)
library(ggraph)
library(tidygraph)

# 適当なネットワークを生成
set.seed(3)  # シード値を固定
G <- sample_gnp(30, 0.1)

g_tidy <- as_tbl_graph(G, directed = FALSE) |>
  activate(nodes) |>
  mutate(
    betweenness = centrality_betweenness(),
    closeness = centrality_closeness(),
    eigen = centrality_eigen(),
    pagerank = centrality_pagerank()
    )

# 描画
p1 <- g_tidy |>
  ggraph() +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(aes(color = betweenness), size = 5) +
  scale_color_viridis(option = "turbo") +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "媒介中心性")

p2 <- g_tidy |>
  ggraph() +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(aes(color = closeness), size = 5) +
  scale_color_viridis(option = "turbo") +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "近接中心性")

p3 <- g_tidy |>
  ggraph() +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(aes(color = eigen), size = 5) +
  scale_color_viridis(option = "turbo") +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "固有値中心性")

p4 <- g_tidy |>
  ggraph() +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(aes(color = pagerank), size = 5) +
  scale_color_viridis(option = "turbo") +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "ページランク")

{p1|p2}/{p3|p4}

7.3.2 様々なネットワーク指標の計算例

7.3.2.1

library(conflicted)
library(tidyverse)
library(igraph)
library(viridis)
library(ggraph)
library(tidygraph)

# ネットワークを生成
n <- 49
m <- 1

dim <- sqrt(n)
G_square <- graph.lattice(length = dim, dim = 2, directed = FALSE) |>
  as_tbl_graph(directed = FALSE)
G_random <- erdos.renyi.game(n, 0.2, directed = FALSE) |>
  as_tbl_graph(directed = FALSE)
G_ba <- barabasi.game(n, m, directed = FALSE) |>
  as_tbl_graph(directed = FALSE)

p1 <- G_square |>
  ggraph(layout = "igraph", algorithm ="grid") +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(color = "#F8766D", size = 5) +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "正方格子")
  
p2 <- G_random |>
  ggraph(layout = "igraph", algorithm ="randomly") +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(color = "#00BA38", size = 5) +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "ランダムネットワーク")

p3 <- G_ba |>
  ggraph(layout = "igraph", algorithm ="fr") +
  geom_edge_link(color = "black", alpha = 0.5) +
  geom_node_point(color = "#619CFF", size = 5) +
  theme_void() +
  theme(aspect.ratio = 1, legend.position = "none") +
  labs(title = "BAネットワーク")

p1 + p2 + p3

7.3.2.2

library(conflicted)
library(tidyverse)
library(igraph)
# library(ggraph)

# ネットワークを生成
n <- 49
m <- 1
dim <- sqrt(n)
G_square <- graph.lattice(length = dim, dim = 2, directed = FALSE)
G_random <- erdos.renyi.game(n, 0.2, directed = FALSE)
G_ba <- barabasi.game(n, m, directed = FALSE)

# ネットワークの指標を計算
compute_network_metrics <- function(G) {
  return(data.frame(
      平均次数 = mean(degree(G)),
      平均最短経路長 = mean_distance(G, directed = FALSE),
      クラスター係数 = transitivity(G, type = "average"),
      同類選択性 = assortativity_degree(G, directed = FALSE)
    ))
}

# 描画
network_names <- c("格子", "ランダム", "BA")

list(G_square, G_random, G_ba) |>
  map(\(x) compute_network_metrics(x)) |>
  list_rbind() |>
  mutate(network_names = factor(network_names, levels = network_names)) |>
  pivot_longer(!network_names) |>
  mutate(name = factor(name, levels = c("平均次数", "平均最短経路長", "クラスター係数", "同類選択性"))) |>
  ggplot(aes(x = network_names, y = value, fill = network_names, label = round(value, 2))) +
  geom_col() +
  geom_text(vjust = -0.5) +
  facet_wrap(vars(name), ncol = 4, scale = "free") +
  theme(legend.position="none", axis.title = element_blank())

第7章はここまで。