第1章

x <- c(5.0, 4.9, 5.3, 5.2, 3.9, 4.5, 6.1, 4.8, 5.2, 4.3)
mu <- 5
sigma <- 3
n <- length(x)
sqrt(n) * (mean(x)-mu) / sigma   ## Zの値
## [1] -0.0843274
qnorm(0.025)                     ## 棄却域との境界(下限)
## [1] -1.959964
qnorm(0.975)                     ## 棄却域との境界(上限)
## [1] 1.959964

第2章

kruskal <- function(w) {     # Prim法
    # 行列の次元取得
    p <- ncol(w)
    # 上三角行列を対称行列に変換
    for (i in 1:(p - 1)) {
        for (j in (i + 1):p) {
            w[j, i] <- w[i, j]
        }
    }
    # 対角要素を 0 に設定
    for (i in 1:p) {
        w[i, i] <- 0
    }
    # 初期化
    parent <- array(dim = p)  # 各ノードの親を記録する配列
    u <- rep(0, p)  # 訪問済みノードを記録 (0: 未訪問, 1: 訪問済み)
    v <- rep(-Inf, p)  # 各ノードの最大重み (最大極大木のため)
    # Uに最初に含まれる頂点jの設定
    j <- 1
    parent[1] <- -1  # 根の親を -1 に設定 (ルート)
    # Prim法による最大極大森の構築
    while (j <= p) {
        u[j] <- 1  # ノード j を訪問済みに設定
        # 隣接ノードの重みを更新
        for (i in 1:p) {
            if (u[i] == 0 && w[j, i] > v[i]) {
                v[i] <- w[j, i]
                parent[i] <- j
            }
        }
        # Uの中にある親iとの重みw[i,j]が最大のV＼Uの要素jを見出す
        max_weight <- 0
        for (i in 1:p) {
            if (u[i] == 0 && v[i] > max_weight) {
                j <- i
                max_weight <- v[i]
            }
        }
        # 極大木にはならない場合、次のjを他の連結部分から見出す
        if (max_weight <= 0) {
            j <- 1
            while (j <= p && u[j] == 1) {
                j <- j + 1
            }
            # すべてu[j]=1であれば、ここでjがpを超えている。
            # 他に連結部分があれば、そのjを新しい根とする。
            if (j <= p) {
                parent[j] <- -1
            }
        }
    }
    # 結果のエッジペアを作成
    pair.1 <- NULL
    pair.2 <- NULL
    for (j in 1:p) {
        if (parent[j] != -1) {
            pair.1 <- c(pair.1, parent[j])
            pair.2 <- c(pair.2, j)
        }
    }
    return(cbind(pair.1, pair.2))  # エッジリストを返す
}

library(igraph)
p <- 10
w <- matrix(rnorm(p*p)**2, p, p)
edgelist <- kruskal(w)
plot(graph_from_edgelist(edgelist), directed=FALSE)

第6章

AIC.1 <- function(k){
  m <- length(k)
  n <- sum(k)
  return(sum(k*log (n/k))+m-1)
}

BIC.1 <- function(k){
  m <- length(k)
  n <- sum(k)
  return(sum(k*log (n/k))+(m-1)/2*log(n))
}

n <- 200
p <- 0.25
x <- rbinom(n, 1, p)
k <- as.vector(table(x))
AIC.1(k)
## [1] 107.6328
BIC.1(k)
## [1] 109.2819
sigma2 <- function(X, y){
  n <- length(y)
  beta <- solve(t(X)%*%X) %*% t(X) %*% y
  return( sum((y-X%*%beta)^2) / n)
}
AIC.2 <- function(X, y){
  n <- length(y)
  p <- ncol(X)
  sigma2 <- sigma2(X, y)
  # AICの公式に代入する
  value <- (n/2)*log(sigma2) + n/2*log(2*pi*exp(1)) + p + 1
  return(value)
}
BIC.2 <- function(X, y){
  n <- length(y)
  p <- ncol(X)
  # 行列演算でbetaの値を得る
  sigma2 <- sigma2(X,y)
  # BICの公式に代入する
  value <- (n/2)*log(sigma2) + n/2*log(2*pi*exp(1)) + (p+1)/2*log(n)
  return(value)
}
AIC.3 <- function(X, y) length(y) * log(sigma2(X, y)) + 2*ncol(X)
BIC.3 <- function(X, y) length(y) * log(sigma2(X, y)) + ncol(X)*log(length(y))

n <- 100
p <- 3
X <- matrix(rnorm(n*p), n, p)
y <- rnorm(n)
mu <- rnorm(p)
Sigma <- diag(p)
AIC.2(X, y)
## [1] 147.9021
BIC.2(X, y)
## [1] 153.1124
AIC.3(X, y)
## [1] 10.01642
BIC.3(X, y)
## [1] 17.83193
library(MASS)
df <- Boston
# 14変数あるうちの離散の説明変数と目的う変数を除外した
X <-as.matrix(df[,c(1,3,5,6,7,8,10,11,12,13)])
y <- df[[14]]
n <- nrow(X)
p <- ncol(X)
X <- cbind(rep(1,n),X) # 切片を含めた
AIC.seq <- NULL
BIC.seq <- NULL;
for(k in 1:p){
  # combn(S,k)で集合Sの中の大きさkの部分集合を列にもつ行列を生成
  T <- combn(1:p,k)
  m <- ncol(T)
  # 各kで、sigma2を最小にする変数の組合わせを選択して、それに対してのAIC/BICを計算。最小値の初期値として S.minに無限大をおいた。
  S.min <- Inf
  for(j in 1:m){
    q <- c(1, T[,j]+1) # 切片(1列目)と選択された変数の列(1列右にずらした)
    S <- sigma2(X[,q],y)
    if(S < S.min) {
      S.min <- S
      q.min <- q
    }
  }
  AIC.seq <- c(AIC.seq, AIC.3(X[,q.min],y))
  BIC.seq <- c(BIC.seq, BIC.3(X[,q.min],y))
}
plot(1:p, ylim=c(min(AIC.seq),max(BIC.seq)), type="n",
     xlab="変数の個数", ylab="AIC/BICの値")
lines(AIC.seq,col="red"); lines(BIC.seq,col="blue")
legend("topright",legend=c("AIC","BIC"), col=c("red","blue"), lwd=1, cex =.8)

Q.1 <- function(k,a){
  m <- length(k)
  log.q <- 0
  for(j in 2:m) {
    log.q <- log.q + log( beta(k[j-1]+a[j-1], sum(k[j:m]+a[j:m])) )
                   - log( beta(a[j-1], sum(a[j:m])))
  }
  return(-log.q)
}

n <- 100
m <- 4
prob <- rep(1/m, m)
k <- as.vector(rmultinom(1, n, prob))
a <- rep(1,m)
Q.1(k,a)
## [1] 144.8654
n <- 100
x <- rbinom(n, 1, 0.5)
y <- rbinom(n, 1, 0.25 )
k <- as.vector(table(x,y))
a <- rep(0.5, 4)
Q.1(k,a)
## [1] 123.5134
# 必要なパッケージを読み込み
library(CholWishart)  # lmvgamma 関数を利用するために必要

# 入力:
#   x       : 観測データ（n行p列の行列形式）
#   kappa.0 : 事前分布のスケールパラメータ
#   nu.0    : 事前分布の自由度
#   mu.0    : 事前分布の平均ベクトル
#   Lambda.0: 事前分布の共分散行列
# 出力:
#   value   : 周辺尤度の対数値

Q.3 <- function(x, kappa.0, nu.0, mu.0, Lambda.0) {
  x <- as.matrix(x)
  n <- nrow(x)
  d <- ncol(x)

  kappa.n <- kappa.0 + n
  nu.n <- nu.0 + n

  x.bar <- colMeans(x)
  S <- t(x - matrix(x.bar, n, d, byrow = TRUE)) %*%
        (x - matrix(x.bar, n, d, byrow = TRUE))
  mu.diff <- matrix(x.bar - mu.0, nrow = d)
  M <- (kappa.0 * n) / (kappa.0 + n) * (mu.diff %*% t(mu.diff))
  Lambda.n <- Lambda.0 + S + M
  value <- n * d / 2 * log(pi) -
    lmvgamma(nu.n / 2, d) + lmvgamma(nu.0 / 2, d) -
    (nu.0 / 2) * log(det(Lambda.0)) + (nu.n / 2) * log(det(Lambda.n)) -
    (d / 2) * log(kappa.0 / kappa.n)
  return(value)
}

library(MASS)
library(CholWishart)
kappa.0 <- 1
nu.0 <- 2
mu.0 <- c(0, 0)
Lambda.0 <- diag(2)
n <- 100
x <- mvrnorm(n, c(0,0), matrix(c(1,0,0,1), 2, 2))
a <- rnorm(2)
A <- a %*% t(a)
y <- mvrnorm(n, c(0,0), A)
Q.3(x, kappa.0, nu.0, mu.0, Lambda.0)
## [1] 266.1074
Q.3(y, kappa.0, nu.0, mu.0, Lambda.0)
## [1] 70.9682
# Bostonデータを準備
df <- Boston
X <- as.matrix(df[, c(1, 3, 5, 6, 7, 8, 10, 11, 12, 13)])  # 説明変数10個
Y <- as.matrix(df[, 14])  # 目的変数
n <- nrow(X)
X <- X[1:n, ]
Y <- Y[1:n, ]
p <- ncol(X)

# スコア計算 (Q3(X,Y) - Q3(X)) を変数数ごとに最小選択
Q3.seq <- numeric(p)
for (k in 1:p) {
  T <- combn(1:p, k)
  m <- ncol(T)
  min_score <- Inf
  for (j in 1:m) {
    idx <- T[, j]
    X.sub <- X[, idx, drop = FALSE]
    mu.0 <- rep(0, ncol(X.sub))
    Lambda.0 <- diag(ncol(X.sub))
    q3_x <- Q.3(X.sub, 1, ncol(X.sub) + 2, mu.0, Lambda.0)
    X.sub.y <- cbind(X.sub, Y)
    mu.0.y <- rep(0, ncol(X.sub.y))
    Lambda.0.y <- diag(ncol(X.sub.y))
    q3_xy <- Q.3(X.sub.y, 1, ncol(X.sub.y) + 2, mu.0.y, Lambda.0.y)
    score <- q3_xy - q3_x
    if (score < min_score) {
      min_score <- score
    }
  }
  Q3.seq[k] <- min_score
}
# プロット
plot(1:p, Q3.seq, type = "b", col = "forestgreen", lwd = 2,
     xlab = "選択された変数の数", ylab = expression(Q[3](X,Y) - Q[3](X)),
     main = "周辺尤度差 Q.3(X,Y)-Q.3(X) の変数選択評価")
grid()

第7章

library(bnlearn)
library(igraph)
# alarmデータの読み込み
df <- alarm
p <- ncol(df)
# ノード名を "1", "2", ..., "p" に置き換える（bnlearn用）
colnames(df) <- as.character(1:p)
# 構造学習（Hill Climbing + BIC）
alarm.hc <- hc(df, score = "bic")
# arcs 関数で (from, to) 辺の情報を取得（文字列型）
mat <- arcs(alarm.hc)
# 文字列を整数に変換（igraphで使うため）
u <- as.integer(mat[,1])  # 親ノード
v <- as.integer(mat[,2])  # 子ノード
edges <- cbind(u, v)      # 辺の行列

X.1 <- c(2, 2, 2, 2, 1, 2, 2, 2, 1, 2)
X.2 <- c(2, 2, 1, 1, 2, 1, 1, 2, 2, 2)
X.3 <- c(1, 1, 2, 2, 3, 1, 3, 1, 1, 1)

tab.1 <- ftable(X.1)
tab.2 <- ftable(X.2)
tab.3 <- ftable(X.3)
tab.1.2 <- ftable(X.1, X.2)
tab.1.3 <- ftable(X.1, X.3)
tab.2.3 <- ftable(X.2, X.3)
tab.1.2.3 <- ftable(X.1, X.2, X.3)

h <- function(tab) {
    m <- nrow(tab)   # 行数（例: 条件変数 X の取りうる値の個数）
    r <- ncol(tab)   # 列数（例: 目的変数 Y の取りうる値の個数）
    h <- 0           # エントロピーの初期化
    for(i in 1:m) {
        ss <- sum(tab[i, ])  # 行 i の合計（= P(X = i) に比例）
        if(ss > 0) {
            for(j in 1:r) {
                if(tab[i, j] > 0) {
                    # 条件付き確率 p(Y = j | X = i) に対する項
                    h <- h - tab[i, j] * log(tab[i, j] / ss)
                }
            }
        }
    }
    return(h)
}

c(h(tab.1), h(tab.2), h(tab.3), h(tab.1.2), h(tab.1.3), h(tab.2.3), h(tab.1.2.3))
## [1] 5.004024 6.730117 9.502705 5.545177 8.588343 6.862250 5.545177
IC <- function(X, y, proc= "AIC"){
    n <- nrow(X)
    tab <- ftable(data.frame(X,y))
    m <- nrow(tab)
    r <- ncol(tab)
    d <- m*(r-1)
    if(proc=="AIC")    IC <- h(tab)+d else
        IC <- h(tab)+d*log(n)/2
    return(IC)
}

IC.min <- function(X, y, proc = "AIC") {
    # X: 説明変数の行列（n × p）
    # y: 応答変数のベクトル（長さ n）
    # proc: "AIC" または "BIC" のどちらを使用するか（文字列）
    p <- ncol(X)      # 説明変数の数
    n <- length(y)    # データの個数（サンプル数）
    L <- 2^p          # 説明変数の全部分集合の数（2^p 通り）
    # 全変数を使わないモデル（定数項のみのモデル）の IC（AIC/BIC の共通部分）
    IC.min <- n * log(sum((y - mean(y))^2) / n)  # 定数項のみのモデルの対数尤度から計算
    set.min <- NULL   # 最小ICを与える変数集合（初期値）
    if (p == 0) {
        # 説明変数がない場合、定数モデルしかないので即座に返す
        return(list(value = IC.min, set = set.min))
    }
    # i は 2 から 2^p まで（i = 1 は空集合に対応するのでスキップ）
    for (i in 2:L) {
        set <- decimal.to.set(i - 1)  # i-1 を2進数に変換し、どの変数を選ぶかを取得
        # 部分集合に対応する説明変数を使ってモデルを適合させ、ICを計算
        if (proc == "AIC") {
            value <- AIC.3(as.matrix(X[, set]), y)
        } else {
            value <- BIC.3(as.matrix(X[, set]), y)
        }
        # ICがこれまでの最小より小さい場合は、値と集合を更新
        if (value < IC.min) {
            IC.min <- value
            set.min <- set
        }
    }
    # 最小ICとそのときの変数集合を返す
    return(list(value = IC.min, set = set.min))
}

library(MASS)
df <- Boston
index <- c(1,3,5,6,7,8,10,11,12,13)
X <- as.matrix(df[, index])
y <- df[[14]]

decimal.to.set <- function(x) if(x==0) NULL else
    c(decimal.to.set(x-2^floor(log2(x))), floor(log2(x))+1)

IC.min(X, y, "AIC")
## $value
## [1] 1655.203
## 
## $set
## [1]  1  3  4  6  8  9 10
p <- ncol(X)
parent <- list()
parent[[1]] <- NULL
for(j in 2:p)parent[[j]] <- IC.min(as.matrix(X[, 1:(j-1)]), X[,j], "BIC")$set
parent
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## [1] 2 3 4
## 
## [[6]]
## [1] 2 4 5
## 
## [[7]]
## [1] 1 2 3 6
## 
## [[8]]
## NULL
## 
## [[9]]
## [1] 1 3 4 6 7 8
## 
## [[10]]
## [1] 1 3 4 5 6 8
# グラフ描画のためのライブラリ
library(igraph)
# 有向辺を格納する行列（最大1000個分を確保しておく）
edge <- matrix(ncol = 2, nrow = 1000)
# 辺のカウント用
k <- 0
p <- 10
# 各変数とその親変数との間に有向辺を張る（j → i）
for (i in 1:p) {
  for (j in parent[[i]]) {
    k <- k + 1
    edge[k, ] <- c(j, i)  # j → i の有向辺を登録
  }
}
# 有向辺のリストから igraph オブジェクトを作成
g <- graph_from_edgelist(edge[1:k, ], directed = TRUE)
# 各ノード（変数）にラベルを付ける
# ラベルは変数名の頭文字（例: crim → "C", indus → "I", ...）
V(g)$label <- c("C", "I", "N", "R", "A", "R", "T", "P", "B", "L")
# グラフの描画（デフォルトスタイル）
layout <- layout_with_fr(g)
plot(g,
     layout = layout,
     vertex.label.cex = 1,    # ラベルの文字サイズ
     vertex.size = 30,          # ノードの円のサイズ（単位: ピクセル程度）
     vertex.label.color = "black", # ラベルの色
     vertex.color = "yellow",     # ノードの円の色
     edge.arrow.size = 0.5          # 矢印の大きさ
)

# データ生成
n <- 100 # 標本数
# xは1/4の確率で1になる二項分布から生成
x <- rbinom(n, 1, 1/4)
k.x <- as.vector(table(x)) # xのカテゴリごとの出現回数
# yは1/2の確率で1になる二項分布から生成
y <- rbinom(n, 1, 1/2)
k.y <- as.vector(table(y)) # yのカテゴリごとの出現回数
# zはxとyの組み合わせ（カテゴリ数増加）
z <- x + 2 * y
k.z <- as.vector(table(z)) # zのカテゴリごとの出現回数
# 独立性の判定
# xとyの独立性が成り立つ場合、結果は負の値になるはず
result <- Q.1(k.x, rep(1/2, 2)) + Q.1(k.y, rep(1/2, 2)) - Q.1(k.z, rep(1/2, 4))
# 結果を出力
print(result)
## [1] -0.2424068
x <- rbinom(n, 1, 1/2)
k.x <- as.vector(table(x))
y <- (rbinom(n, 1, 1/3) + x) %% 2
k.y <- as.vector(table(y))
z <- x + 2 * y
k.z <- as.vector(table(z))

# データの生成
n <- 100  # サンプル数
a <- rnorm(2) * 10              # ランダムな係数
A <- a %*% t(a)                 # 共分散行列の生成
z <- mvrnorm(n, c(0, 0), A)     # 多変量正規分布データの生成
# x, yに分解
x <- z[, 1]  # zの1列目
y <- z[, 2]  # zの2列目
# 関数の実行例
Q.3(x, 2, 2, c(0), matrix(1, 1, 1))
## [1] 376.8391
Q.3(y, 2, 2, c(0), matrix(1, 1, 1))
## [1] 361.0308
Q.3(z, 2, 2, c(0, 0), diag(2))
## [1] 318.1479
Q.1(k.x, rep(1/2, 2)) + Q.1(k.y, rep(1/2, 2)) - Q.1(k.z, rep(1/2, 4))
## [1] 1.320075
z <- mvrnorm(n, c(0,0), diag(2))

Q.3(x, 2, 2, c(0), matrix(1, 1, 1)) +
Q.3(y, 2, 2, c(0), matrix(1, 1, 1)) -
Q.3(z, 2, 2, c(0, 0), diag(2))
## [1] 431.1244
I.n <- function(x, y) {
    z <- table(x, y)  # X, Y の同時度数表
    size <- dim(z)  # 表のサイズを取得
    n <- length(x)  # サンプルサイズ
    S <- 0  # 相互情報量の初期値
    for (i in 1:size[1]) {
        u <- sum(z[i,])  # X = i の頻度
        for (j in 1:size[2]) {
            v <- sum(z[, j])  # Y = j の頻度
            if (z[i, j] != 0) {
                S <- S + z[i, j] / n * log( z[i, j] * n / (u * v) )
                # 相互情報量の計算
            }
        }
    }
    return(S)
}

library(igraph)
library(bnlearn)
# データの取得
df <- asia
n <- nrow(df)  # サンプル数
p <- ncol(df)  # 変数の個数
# データを行列形式に変換
x <- matrix(nrow = n, ncol = p)
for (i in 1:p) {
    x[, i] <- df[[i]]
}
# 相互情報量行列の作成
w <- matrix(0, nrow = p, ncol = p)  # 初期化
for (i in 1:(p - 1)) {
    for (j in (i + 1):p) {
        w[i, j] <- I.n(x[, i], x[, j])  # 相互情報量を計算
        w[j, i] <- w[i, j]
    }
}
# Chow-Liu 木の構築（最大重み全域木）
edge.list <- kruskal(w)  # Kruskalアルゴリズム
print(edge.list)  # 結果のエッジリストを表示
##      pair.1 pair.2
## [1,]      4      2
## [2,]      1      3
## [3,]      6      4
## [4,]      2      5
## [5,]      3      6
## [6,]      6      7
## [7,]      5      8
# グラフを描画
g <- graph_from_edgelist(edge.list, directed = FALSE)  # igraph を使用
V(g)$label <- colnames(df)  # ノードにラベルを付与
plot(g, main = "Chow-Liu アルゴリズム(I.n)")  # グラフを描画

J.n <- function(x,y){
  x <- as.numeric(x)
  alpha <- length(table(x))
  y <- as.numeric(y)
  beta <- length(table(y))
  z <- (y-1) * alpha + x
  value <- Q.4(x, alpha) + Q.4(y, beta) - Q.4(z, alpha*beta)
  return(max(value,0))
}

Q.4 <- function(x, alpha){
  x <- as.numeric(x)
  n <- length(x)
  nc <- 0
  cc <- rep(0, alpha)
  log.q <- 0
  for(i in 1:n){
    log.q <- log.q - log ( (cc[x[i]]+0.5) / (nc + alpha/2))
    cc[x[i]] <- cc[x[i]] + 1
    nc <- nc +1
  }
  return(log.q / n)
}

J.n <- function(x,y){
  n <- length (x)
  size <- dim(table(x,y))
  value <- I.n(x,y) - (size[1]-1) * (size[2]-1) / 2 / n * log(n)
  return(max(value, 0))
}

I.n <- function(x,y) {
  rho <- cov(x, y) / sqrt(var(x) * var(y))
  return(-0.5 * log(1-rho^2))
}

J.n <- function(x,y){
  n <- length(x)
  value <- I.n(x,y) - 0.5 * log (n) / n
  return (max(value, 0))
}