5.2 グラフィカルLasso
inner.prod <- function(x, y) {
return(sum(x * y))
} ## 既出
soft.th <- function(lambda, x) {
return(sign(x) * pmax(abs(x) - lambda, 0))
} ## 既出
graph.lasso <- function(s, lambda = 0) {
W <- s
p <- ncol(s)
beta <- matrix(0, nrow = p - 1, ncol = p)
beta.out <- beta
eps.out <- 1
while (eps.out > 0.01) {
for (j in 1:p) {
a <- W[-j, -j]
b <- s[-j, j]
beta.in <- beta[, j]
eps.in <- 1
while (eps.in > 0.01) {
for (h in 1:(p - 1)) {
cc <- b[h] - inner.prod(a[h, -h], beta[-h, j])
beta[h, j] <- soft.th(lambda, cc) / a[h, h]
}
eps.in <- max(beta[, j] - beta.in)
beta.in <- beta[, j]
}
W[-j, j] <- W[-j, -j] %*% beta[, j]
}
eps.out <- max(beta - beta.out)
beta.out <- beta
}
theta <- matrix(nrow = p, ncol = p)
for (j in 1:p) {
theta[j, j] <- 1 / (W[j, j] - W[j, -j] %*% beta[, j])
theta[-j, j] <- -beta[, j] * theta[j, j]
}
return(theta)
}
例47
library(MultiRNG)
## Warning: package 'MultiRNG' was built under R version 4.0.3
Theta <- matrix(c( 2, 0.6, 0, 0, 0.5, 0.6, 2, -0.4, 0.3, 0,
0, -0.4, 2, -0.2, 0, 0, 0.3, -0.2, 2, -0.2,
0.5, 0, 0, -0.2, 2),
nrow = 5)
Sigma <- solve(Theta)
meanvec <- rep(0, 5)
dat <- draw.d.variate.normal(no.row = 20, d = 5,
mean.vec = meanvec, cov.mat = Sigma)
# 平均mean.vec，共分散行列cov.mat，サンプル数no.row，変数の個数dから
# サンプル行列を生成
s <- t(dat) %*% dat / nrow(dat)
Theta
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.0 0.6 0.0 0.0 0.5
## [2,] 0.6 2.0 -0.4 0.3 0.0
## [3,] 0.0 -0.4 2.0 -0.2 0.0
## [4,] 0.0 0.3 -0.2 2.0 -0.2
## [5,] 0.5 0.0 0.0 -0.2 2.0
graph.lasso(s)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.5026572 0.9818782 0.60009124 -2.0466927 1.33991063
## [2,] 0.9751595 1.9353862 -0.23683861 0.1526158 0.10016728
## [3,] 0.5963550 -0.2364995 2.00318046 0.2183919 0.09857318
## [4,] -2.0529957 0.1479875 0.21616558 3.4705123 -0.92298628
## [5,] 1.3409907 0.1029790 0.09968043 -0.9214949 1.87734886
graph.lasso(s, lambda = 0.015)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.9201795 0.83909586 0.4968786 -1.6663040 1.00974904
## [2,] 0.8419598 1.86659378 -0.2258595 0.1875745 0.01679017
## [3,] 0.4997081 -0.22613987 1.9474979 0.2293072 0.02867895
## [4,] -1.6599229 0.18933948 0.2312845 3.1619637 -0.69464559
## [5,] 1.0079087 0.01586028 0.0280105 -0.6952598 1.73576457
graph.lasso(s, lambda = 0.03)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.5384878 0.7691536 0.4401780 -1.3900787 0.7967628
## [2,] 0.7678327 1.8114040 -0.2103822 0.1852676 0.0000000
## [3,] 0.4410263 -0.2110721 1.9023682 0.2157557 0.0000000
## [4,] -1.3890172 0.1847097 0.2163867 2.9235424 -0.5134782
## [5,] 0.7972078 0.0000000 0.0000000 -0.5137533 1.6430676
graph.lasso(s, lambda = 0.05)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.1732901 0.6937337 0.3898361 -1.1235715 0.5997789
## [2,] 0.6921656 1.7451215 -0.1881856 0.1691241 0.0000000
## [3,] 0.3893990 -0.1888305 1.8497121 0.1808854 0.0000000
## [4,] -1.1237883 0.1683662 0.1809569 2.6903188 -0.3216729
## [5,] 0.6006711 0.0000000 0.0000000 -0.3218116 1.5604271
例48
library(glasso)
## Warning: package 'glasso' was built under R version 4.0.3
solve(s)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.4996092 0.9759424 0.59616546 -2.0504223 1.33984091
## [2,] 0.9759424 1.9340685 -0.23601615 0.1513602 0.10110823
## [3,] 0.5961655 -0.2360162 2.00218014 0.2187106 0.09831414
## [4,] -2.0504223 0.1513602 0.21871064 3.4723473 -0.92260766
## [5,] 1.3398409 0.1011082 0.09831414 -0.9226077 1.87722288
glasso(s, rho = 0)
## Warning in glasso(s, rho = 0): With rho=0, there may be convergence problems if the input matrix
## is not of full rank
## $w
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4840149 -0.2783125 -0.1957499 0.25901256 -0.19291794
## [2,] -0.2783125 0.6915990 0.1823482 -0.19050703 0.05820960
## [3,] -0.1957499 0.1823482 0.5950644 -0.15503142 0.02253190
## [4,] 0.2590126 -0.1905070 -0.1550314 0.47706857 0.06798044
## [5,] -0.1929179 0.0582096 0.0225319 0.06798044 0.69949020
##
## $wi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.4996218 0.9759768 0.59618930 -2.0504188 1.33983912
## [2,] 0.9759458 1.9340837 -0.23601549 0.1513473 0.10111612
## [3,] 0.5961611 -0.2360212 2.00218797 0.2187012 0.09831833
## [4,] -2.0504304 0.1513429 0.21869822 3.4723342 -0.92260025
## [5,] 1.3398450 0.1011183 0.09832162 -0.9226041 1.87721935
##
## $loglik
## [1] NA
##
## $errflag
## [1] 0
##
## $approx
## [1] FALSE
##
## $del
## [1] 1.422557e-06
##
## $niter
## [1] 1
glasso(s, rho = 0.015)
## $w
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4990149 -0.26331132 -0.18074581 0.24401373 -0.17791758
## [2,] -0.2633113 0.70659903 0.16734523 -0.17550769 0.07320939
## [3,] -0.1807458 0.16734523 0.61006437 -0.14003150 0.03753097
## [4,] 0.2440137 -0.17550769 -0.14003150 0.49206857 0.05298044
## [5,] -0.1779176 0.07320939 0.03753097 0.05298044 0.71449020
##
## $wi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.6395267 0.780865279 0.46547486 -1.4920883 0.912469153
## [2,] 0.7808639 1.805685316 -0.21995497 0.1935040 0.006634033
## [3,] 0.4654826 -0.219960606 1.88779032 0.2255145 0.022562388
## [4,] -1.4920810 0.193502603 0.22551594 2.9724949 -0.623636388
## [5,] 0.9124746 0.006634729 0.02256159 -0.6236402 1.671194883
##
## $loglik
## [1] NA
##
## $errflag
## [1] 0
##
## $approx
## [1] FALSE
##
## $del
## [1] 1.391914e-06
##
## $niter
## [1] 2
glasso(s, rho = 0.030)
## $w
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.5140149 -0.24830961 -0.16574571 0.22901326 -0.16291838
## [2,] -0.2483096 0.72159903 0.15234577 -0.16050791 0.06508705
## [3,] -0.1657457 0.15234577 0.62506437 -0.12503155 0.03856086
## [4,] 0.2290133 -0.16050791 -0.12503155 0.50706857 0.03798044
## [5,] -0.1629184 0.06508705 0.03856086 0.03798044 0.72949020
##
## $wi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.1262203 0.6756426 0.3920615 -1.1521078 0.6771599
## [2,] 0.6756458 1.7012520 -0.1986007 0.1843966 0.0000000
## [3,] 0.3920730 -0.1986037 1.7926239 0.2020837 0.0000000
## [4,] -1.1521061 0.1843949 0.2020837 2.6322279 -0.4214819
## [5,] 0.6771633 0.0000000 0.0000000 -0.4214826 1.5439959
##
## $loglik
## [1] NA
##
## $errflag
## [1] 0
##
## $approx
## [1] FALSE
##
## $del
## [1] 2.075307e-06
##
## $niter
## [1] 2
glasso(s, rho = 0.045)
## $w
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.5290149 -0.23330907 -0.15074659 0.21401278 -0.14791839
## [2,] -0.2333091 0.73659903 0.13734602 -0.14550749 0.05551342
## [3,] -0.1507466 0.13734602 0.64006437 -0.11003147 0.03282273
## [4,] 0.2140128 -0.14550749 -0.11003147 0.52206857 0.02298044
## [5,] -0.1479184 0.05551342 0.03282273 0.02298044 0.74449020
##
## $wi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.7652114 0.5917274 0.3395463 -0.9198957 0.5187054
## [2,] 0.5917283 1.6115062 -0.1773502 0.1692015 0.0000000
## [3,] 0.3395498 -0.1773516 1.7098944 0.1717574 0.0000000
## [4,] -0.9198980 0.1692008 0.1717572 2.3880891 -0.2766716
## [5,] 0.5187076 0.0000000 0.0000000 -0.2766721 1.4547997
##
## $loglik
## [1] NA
##
## $errflag
## [1] 0
##
## $approx
## [1] FALSE
##
## $del
## [1] 1.202403e-06
##
## $niter
## [1] 2
library(igraph)
## Warning: package 'igraph' was built under R version 4.0.3
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
adj <- function(mat) {
p <- ncol(mat)
ad <- matrix(0, nrow = p, ncol = p)
for (i in 1:(p - 1)) {
for (j in (i + 1):p) {
if (mat[i, j] == 0) {
ad[i, j] <- 0
} else {
ad[i, j] <- 1
}
}
}
g <- graph.adjacency(ad, mode = "undirected")
plot(g)
}
例49
library(glasso)
library(igraph)
df <- read.csv("breastcancer.csv")
w <- matrix(nrow = 250, ncol = 1000)
for (i in 1:1000)
w[, i] <- as.numeric(df[[i]])
x <- w
s <- t(x) %*% x / 250
fit <- glasso(s, rho = 0.75)
sum(fit$wi == 0)
## [1] 992822
y <- NULL
z <- NULL
for (i in 1:999) {
for (j in (i + 1):1000) {
if (fit$wi[i, j] != 0) {
y <- c(y, i)
z <- c(z, j)
}
}
}
edges <- cbind(y, z)
write.csv(edges,"edges.csv")
5.3 疑似尤度を用いたグラフィカルモデルの推定
例50
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.0.3
## Loading required package: Matrix
## Loaded glmnet 4.0-2
df <- read.csv("breastcancer.csv")
n <- 250
p <- 50
w <- matrix(nrow = n, ncol = p)
for (i in 1:p)
w[, i] <- as.numeric(df[[i]])
x <- w[, 1:p]
fm <- rep("gaussian", p)
lambda <- 0.1
fit <- list()
for (j in 1:p)
fit[[j]] <- glmnet(x[, -j], x[, j], family = fm[j], lambda = lambda)
ad <- matrix(0, p, p)
for (i in 1:p) {
for (j in 1:(p - 1)) {
k <- j
if (j >= i)
k <- j + 1
if (fit[[i]]$beta[j] != 0) {
ad[i, k] <- 1
} else {
ad[i, k] <- 0
}
}
}
## ANDの場合
for (i in 1:(p - 1)) {
for (j in (i + 1):p) {
if (ad[i, j] != ad[i, j]) {
ad[i, j] <- 0
ad[j, i] <- 0
}
}
}
u <- NULL
v <- NULL
for (i in 1:(p - 1)) {
for (j in (i + 1):p) {
if (ad[i, j] == 1) {
u <- c(u, i)
v <- c(v, j)
}
}
}
u
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3
## [26] 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6
## [51] 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9
## [76] 9 9 9 9 9 9 10 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12
## [101] 12 12 13 13 13 13 13 13 13 13 13 14 14 15 15 15 16 16 16 16 16 16 16 16 16
## [126] 16 16 17 17 17 17 17 18 18 18 18 18 18 18 18 19 19 19 19 20 20 20 21 21 21
## [151] 21 21 22 22 22 22 22 22 23 23 23 23 23 23 23 24 24 24 24 24 24 25 25 25 25
## [176] 26 26 26 26 26 27 27 27 27 27 28 28 28 28 29 29 29 29 29 30 30 30 30 30 30
## [201] 31 31 31 31 31 31 31 32 32 33 33 34 34 34 35 35 36 36 36 37 37 37 38 38 38
## [226] 38 38 40 40 42 42 43 43 43 44 44 45 45 47 47 47 49
v
## [1] 2 12 13 18 22 23 24 28 41 46 48 50 4 5 10 16 18 21 23 26 27 34 37 43 9
## [26] 10 15 19 20 21 24 29 38 50 5 8 17 18 21 23 46 9 20 21 28 30 32 10 11 15
## [51] 21 24 25 31 33 39 41 45 8 12 14 18 30 36 44 47 12 13 15 17 41 49 10 14 18
## [76] 19 21 29 45 46 50 18 25 28 32 34 37 42 12 16 20 28 42 48 17 18 29 33 35 41
## [101] 43 47 20 22 25 28 39 40 43 45 50 15 16 17 25 45 18 23 26 27 30 31 32 41 46
## [126] 47 49 20 25 29 31 36 21 24 27 35 41 43 48 49 20 25 42 44 27 44 48 23 24 34
## [151] 37 42 23 29 30 35 36 50 25 28 32 35 37 47 50 28 33 35 40 41 42 29 39 45 48
## [176] 39 40 42 46 50 37 40 42 43 49 29 31 47 50 30 39 45 48 50 31 32 45 47 48 50
## [201] 34 36 39 44 47 49 50 33 45 37 49 35 36 45 41 42 38 39 45 38 42 46 39 45 46
## [226] 49 50 41 42 49 50 47 49 50 48 50 46 50 48 49 50 50
adj(ad)
## ORの場合
for (i in 1:(p - 1)) {
for (j in (i + 1):p) {
if (ad[i, j] != ad[i, j]) {
ad[i, j] <- 1
ad[j, i] <- 1
}
}
}
adj(ad)
例51
library(glmnet)
df <- read.csv("breastcancer.csv")
w <- matrix(nrow = 250, ncol = 1000)
for (i in 1:1000)
w[, i] <- as.numeric(df[[i]])
w <- (sign(w) + 1) / 2 ## 2値データに変換する
p <- 50
x <- w[, 1:p]
fm <- rep("binomial", p)
lambda <- 0.15
fit <- list()
for (j in 1:p)
fit[[j]] <- glmnet(x[, -j], x[, j], family = fm[j], lambda = lambda)
ad <- matrix(0, nrow = p, ncol = p)
for (i in 1:p) {
for (j in 1:(p - 1)) {
k <- j
if (j >= i)
k <- j + 1
if (fit[[i]]$beta[j] != 0) {
ad[i, k] <- 1
} else {
ad[i, k] <- 0
}
}
}
for (i in 1:(p - 1)) {
for (j in (i + 1):p) {
if (ad[i, j] != ad[i, j]) {
ad[i, j] <- 0
ad[j, i] <- 0
}
}
}
sum(ad)
## [1] 203
adj(ad)
5.4 JointグラフィカルLasso
# 大きさが3以上でないとgenlassoは稼働しない
b.fused <- function(y, lambda) {
if (y[1] > y[2] + 2 * lambda) {
a <- y[1] - lambda
b <- y[2] + lambda
} else if (y[1] < y[2] - 2 * lambda) {
a <- y[1] + lambda
b <- y[2] - lambda
} else {
a <- (y[1] + y[2]) / 2
b <- a
}
return(c(a, b))
}
# 隣接項だけではなく，離接するすべての値と比較するFused Lasso
fused <- function(y, lambda.1, lambda.2) {
K <- length(y)
if (K == 1) {
theta <- y
} else if (K == 2) {
theta <- b.fused(y, lambda.2)
} else {
L <- K * (K - 1) / 2
D <- matrix(0, nrow = L, ncol = K)
k <- 0
for (i in 1:(K - 1)) {
for (j in (i + 1):K) {
k <- k + 1
D[k, i] <- 1
D[k, j] <- -1
}
}
out <- genlasso(y, D = D)
theta <- coef(out, lambda = lambda.2)
}
theta <- soft.th(lambda.1, theta)
return(theta)
}
# Joint Graphical Lasso
jgl <- function(X, lambda.1, lambda.2) { # Xはリストで与える
K <- length(X)
p <- ncol(X[[1]])
n <- array(dim = K)
S <- list()
for (k in 1:K) {
n[k] <- nrow(X[[k]])
S[[k]] <- t(X[[k]]) %*% X[[k]] / n[k]
}
rho <- 1
lambda.1 <- lambda.1 / rho
lambda.2 <- lambda.2 / rho
Theta <- list()
for (k in 1:K)
Theta[[k]] <- diag(p)
Theta.old <- list()
for (k in 1:K)
Theta.old[[k]] <- diag(rnorm(p))
U <- list()
for (k in 1:K)
U[[k]] <- matrix(0, nrow = p, ncol = p)
Z <- list()
for (k in 1:K)
Z[[k]] <- matrix(0, nrow = p, ncol = p)
epsilon <- 0
epsilon.old <- 1
while (abs((epsilon - epsilon.old) / epsilon.old) > 0.0001) {
Theta.old <- Theta
epsilon.old <- epsilon
## (a)に関する更新
for (k in 1:K) {
mat <- S[[k]] - rho * Z[[k]] / n[k] + rho * U[[k]] / n[k]
svd.mat <- svd(mat)
V <- svd.mat$v
D <- svd.mat$d
DD <- n[k] / (2 * rho) * (-D + sqrt(D ^ 2 + 4 * rho / n[k]))
Theta[[k]] <- V %*% diag(DD) %*% t(V)
}
## (b)に関する更新
for (i in 1:p) {
for (j in 1:p) {
A <- NULL
for (k in 1:K)
A <- c(A, Theta[[k]][i, j] + U[[k]][i, j])
if (i == j) {
B <- fused(A, 0, lambda.2)
} else {
B <- fused(A, lambda.1, lambda.2)
}
for (k in 1:K)
Z[[k]][i, j] <- B[k]
}
}
## (c)に関する更新
for (k in 1:K)
U[[k]] <- U[[k]] + Theta[[k]] - Z[[k]]
## 収束したかどうかの検査
epsilon <- 0
for (k in 1:K) {
epsilon.new <- max(abs(Theta[[k]] - Theta.old[[k]]))
if (epsilon.new > epsilon)
epsilon <- epsilon.new
}
}
return(Z)
}
## (b)の更新の箇所を以下で置き換える。下記単独では動作しない。
for (i in 1:p) {
for (j in 1:p) {
A <- NULL
for (k in 1:K)
A <- c(A, Theta[[k]][i, j] + U[[k]][i, j])
if (i == j) {
B <- A
} else {
B <- soft.th(lambda.1 / rho,A) *
max(1 - lambda.2 / rho / sqrt(norm(soft.th(lambda.1 / rho, A), "2") ^ 2), 0)
}
for (k in 1:K)
Z[[k]][i, j] <- B[k]
}
}
例52
## データ生成と実行
p <- 10
K <- 2
N <- 100
n <- array(dim = K)
for (k in 1:K)
n[k] <- N / K
X <- list()
X[[1]] <- matrix(rnorm(n[k] * p), ncol = p)
for (k in 2:K)
X[[k]] <- X[[k - 1]] + matrix(rnorm(n[k] * p) * 0.1, ncol = p)
## lambda.1,lambda.2を変えて実行してみた
Theta <- jgl(X, 3, 0.01)
par(mfrow = c(1, 2))
adj(Theta[[1]])
adj(Theta[[2]])
