第2章

20


f <- function(x)exp(beta.0+beta*x)/(1+exp(beta.0+beta*x))
beta.0 <- 0
beta.seq <- c(0,0.2,0.5,1,2,10)
m=length(beta.seq)
beta <- beta.seq[1]
plot(f,xlim=c(-10,10),ylim=c(0,1),xlab="x",ylab="y",col=1,main="ロジスティック曲線")
for(i in 2:m){
  beta <- beta.seq[i]; par(new=TRUE)
  plot(f,xlim=c(-10,10),ylim=c(0,1),xlab="",ylab="",axes=FALSE,col=i)
}
legend("topleft",legend=beta.seq,col=1:length(beta.seq),lwd=2,cex=.8)

par(new=FALSE)

22


f <- function(x) x^2-1;
f.  < function(x) 2*x;
curve(f(x),-1,2.5); abline(h = 0, col = "blue")
x <- 2
for(i in 1:100){
  X <- x
  Y <- f(x)
  x <- x - f(x) / f.(x)
  y <- f(x)
  segments(X,Y,x,0); segments(X,Y,X,0, lty=3)
  points(x,0,col="red",pch = 16)
}

## [1] 1


f <- function(z) z[1]^2 + z[2]^2 - 1
f.x <- function(z) 2*z[1]
f.y <- function(z)2* z[2]
g <- function(z) z[1] + z[2]
g.x <- function(z)1
g.y <- function(z) 1
z <- c(1,2)
curve(-x,xlim=c(-4,1),ylim=c(-1,4)); curve(sqrt(1-x^2),xlim=c(-1,1),add = TRUE); curve(-sqrt(1-x^2),xlim=c(-1,1),add = TRUE)
points(z[1],z[2],col="red",pch = 16)
for(i in 1:100){
  Z <- z
  z <- z - solve(matrix(c(f.x(z),f.y(z),g.x(z),g.y(z)),ncol=2,byrow=TRUE))%*%c(f(z),g(z))
  segments(Z[1],Z[2],z[1],z[2],lty = 3); 
  points(z[1],z[2],col="red",pch = 16)
}

##            [,1]
## [1,] -0.7071068
## [2,]  0.7071068

24


N <- 1000
p <- 2
X <- matrix(rnorm(N*p),ncol=p)
X <- cbind(rep(1,N),X)
beta <- rnorm(p+1)
y <- array(N)
s <- as.vector(X%*%beta)
prob <- 1/(1+exp(s))
for(i in 1:N){
  if(runif(1)>prob[i]){
    y[i] <- 1 
  }else {
    y[i] <- -1
  }
}
beta


## [1]  0.5585961 -0.2603251  0.3784372

##  最尤推定  ##
beta <- Inf
gamma <- rnorm(p+1)
while(sum((beta-gamma)^2)>0.001){
  beta <- gamma
  s <- as.vector(X%*%beta)
  v <- exp(-s*y)
  u <- y*v/(1+v)
  w <- v/(1+v)^2
  W <- diag(w)
  z <- s+u/w
  gamma <- as.vector(solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z)
  print(gamma)
}

## [1]  0.5655375 -0.2002013  0.3233685
## [1]  0.6037300 -0.2234176  0.3591788
## [1]  0.6044193 -0.2238991  0.3598999

26


n <- 100
x <- c(rnorm(n)+1, rnorm(n)-1)
y <- c(rep(1,n),rep(-1,n))
train <- sample(1:(2*n),n,replace=FALSE)
df <- data.frame(x,y)
x <- as.matrix(df[train,1])
y <- as.vector(df[train,2])
p <- 1
X <- cbind(1,x)
beta <- 0
gamma <- rnorm(p+1)
while(sum((beta-gamma)^2)>0.001){
  beta <- gamma
  s <- as.vector(X%*%beta)
  v <- exp(-s*y)
  u <- y*v/(1+v);
  w <- v/(1+v)^2; W=diag(w)
  z <- s+u/w
  gamma <- as.vector(solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z)
  print(gamma)
}

## [1] -0.1684353  1.3183688
## [1] -0.1626278  1.3293458


x <- as.matrix(df[-train,1])
y <- as.vector(df[-train,2])
z <- 2*as.integer(beta[1]+x*beta[2]>0)-1
table(y,z)

##     z
## y    -1  1
##   -1 34  9
##   1  11 46

28

変更前

# データ生成
mu.1 <- 2
mu.2 <- 2
sigma.1 <- 2
sigma.2 <- 2
rho.1 <- 0

mu.3 <- -3
mu.4 <- -3
sigma.3 <- 1
sigma.4 <- 1
rho.2 <- -0.8

n <- 100

u <- rnorm(n)
v <- rnorm(n)
x.1 <- sigma.1 * u + mu.1
y.1 <- (rho.1 * u + sqrt(1 - rho.1 ^ 2) * v) * sigma.2 + mu.2

u <- rnorm(n)
v <- rnorm(n)
x.2 <- sigma.3 * u + mu.3
y.2 <- (rho.2 * u + sqrt(1 - rho.2 ^ 2) * v) * sigma.4 + mu.4

# 分布を推定して，境界線を引く
f <- function(w, mu, inv, de) {
  return(-0.5 * t(w - mu) %*% inv %*% (w - mu) - 0.5 * log(de))
}

mu.1 <- c(mean(x.1), mean(y.1))
mu.2 <- c(mean(x.2), mean(y.2))

df <- data.frame(x.1, y.1)  #
mat <- cov(df)              #
inv.1 <- solve(mat)         #
de.1 <- det(mat)            #

df <- data.frame(x.2, y.2)  #
mat <- cov(df)   
inv.2 <- solve(mat)         #
de.2 <- det(mat)            #

f.1 <- function(u, v) {
  return(f(c(u, v), mu.1, inv.1, de.1))
}

f.2 <- function(u, v) {
  return(f(c(u, v), mu.2, inv.2, de.2))
}

u <- v <- seq(-6, 6, length = 50)
n <- length(u)
w <- array(dim = c(n, n))
pi.1 <- 0.5
pi.2 <- 0.5
for (i in 1:n)
  for (j in 1:n)
    w[i, j] <- pi.1 * f.1(u[i], v[j]) - pi.2 * f.2(u[i], v[j])
contour(u, v, w, level = 1)
points(x.1, y.1, col = "red")
points(x.2, y.2, col = "blue")

変更後

# データ生成
mu.1 <- 2
mu.2 <- 2
sigma.1 <- 2
sigma.2 <- 2
rho.1 <- 0

mu.3 <- -3
mu.4 <- -3
sigma.3 <- 1
sigma.4 <- 1
rho.2 <- -0.8

n <- 100

u <- rnorm(n)
v <- rnorm(n)
x.1 <- sigma.1 * u + mu.1
y.1 <- (rho.1 * u + sqrt(1 - rho.1 ^ 2) * v) * sigma.2 + mu.2

u <- rnorm(n)
v <- rnorm(n)
x.2 <- sigma.3 * u + mu.3
y.2 <- (rho.2 * u + sqrt(1 - rho.2 ^ 2) * v) * sigma.4 + mu.4

# 分布を推定して，境界線を引く
f <- function(w, mu, inv, de) {
  return(-0.5 * t(w - mu) %*% inv %*% (w - mu) - 0.5 * log(de))
}

mu.1 <- c(mean(x.1), mean(y.1))
mu.2 <- c(mean(x.2), mean(y.2))

df <- data.frame(x.1, y.1)  #
mat <- cov(df)              #
inv.1 <- solve(mat)         #
de.1 <- det(mat)            #

inv.2 <- solve(mat)         #
de.2 <- det(mat)            #

f.1 <- function(u, v) {
  return(f(c(u, v), mu.1, inv.1, de.1))
}

f.2 <- function(u, v) {
  return(f(c(u, v), mu.2, inv.2, de.2))
}

u <- v <- seq(-6, 6, length = 50)
n <- length(u)
w <- array(dim = c(n, n))
pi.1 <- 0.5
pi.2 <- 0.5
for (i in 1:n)
  for (j in 1:n)
    w[i, j] <- pi.1 * f.1(u[i], v[j]) - pi.2 * f.2(u[i], v[j])
contour(u, v, w, level = 1)
points(x.1, y.1, col = "red")
points(x.2, y.2, col = "blue")

29

f <- function(w, mu, inv, de) {
  return(-0.5 * (w - mu) %*% inv %*% t(w - mu) - 0.5 * log(de))
}

df <- iris
df[[5]] <- c(rep(1, 50), rep(2, 50), rep(3, 50))
n <- nrow(df)
train <- sample(1:n, n / 2, replace = FALSE)
test <- setdiff(1:n, train)
mat <- as.matrix(df[train, ])
mu <- list()
covv <- list()
for (j in 1:3) {
  x <- mat[mat[, 5] == j, 1:4]
  mu[[j]] <- c(mean(x[, 1]), mean(x[, 2]), mean(x[, 3]), mean(x[, 4]))
  covv[[j]] <- cov(x)
}

g <- function(v, j) {
  return(f(v, mu[[j]], solve(covv[[j]]), det(covv[[j]])))
}

z <- array(dim = length(test))
for (i in test) {
  u <- as.matrix(df[i, 1:4])
  a <- g(u, 1)
  b <- g(u, 2)
  c <- g(u, 3)
  if (a < b) {
    if (b < c) {
      z[i] <- 3
    } else {
      z[i] <- 2
    }
  } else {
    if (a < c) {
      z[i] <- 3
    } else {
      z[i] <- 1
    }
  }
}
table(z[test], df[test, 5])

##    
##      1  2  3
##   1 27  0  0
##   2  0 20  0
##   3  0  2 26

f <- function(w, mu, inv, de) {
  return(-0.5 * (w - mu) %*% inv %*% t(w - mu) - 0.5 * log(de))
}

df <- iris[-76:-125,]
df[[5]] <- c(rep(1, 50), rep(2, 25), rep(3, 25))
n <- nrow(df)
train <- sample(1:n, n / 2, replace = FALSE)
test <- setdiff(1:n, train)
mat <- as.matrix(df[train, ])
mu <- list()
covv <- list()
for (j in 1:3) {
  x <- mat[mat[, 5] == j, 1:4]
  mu[[j]] <- c(mean(x[, 1]), mean(x[, 2]), mean(x[, 3]), mean(x[, 4]))
  covv[[j]] <- cov(x)
}

g <- function(v, j) {
  return(f(v, mu[[j]], solve(covv[[j]]), det(covv[[j]])))
}

z <- array(dim = length(test))
for (i in test) {
  u <- as.matrix(df[i, 1:4])
  a <- g(u, 1) + log(0.5)
  b <- g(u, 2) + log(0.25)
  c <- g(u, 3) + log(0.25)
  if (a < b) {
    if (b < c) {
      z[i] <- 3
    } else {
      z[i] <- 2
    }
  } else {
    if (a < c) {
      z[i] <- 3
    } else {
      z[i] <- 1
    }
  }
}
table(z[test], df[test, 5])

##    
##      1  2  3
##   1 22  0  0
##   2  0 13  0
##   3  0  2 13

30

knn.1 <- function(x, y, z, k) {
  x <- as.matrix(x)
  n <- nrow(x)
  p <- ncol(x)
  dis <- array(dim = n)
  for (i in 1:n)
    dis[i] <- norm(z - x[i, ], "2")
  S <- order(dis)[1:k]  # 距離の小さいk個のiの集合
  u <- sort(table(y[S]), decreasing = TRUE)  # 最頻のy[i]と頻度
  v <- names(u)[1]  # 最頻のy[i]
  return(v)
}

knn <- function(x, y, z, k) {
  n <- nrow(z)
  w <- array(dim = n)
  for (i in 1:n)
    w[i] <- knn.1(x, y, z[i, ], k)
  return(w)
}

df <- iris
n <- 150
train <- sample(1:n, n / 2, replace = FALSE)
test <- setdiff(1:n, train)
x <- as.matrix(df[train, 1:4])
y <- as.vector(df[train, 5])
z <- as.matrix(df[test, 1:4])
ans <- as.vector(df[test, 5])
w <- knn(x, y, z, k = 3)
table(w, ans)

##             ans
## w            setosa versicolor virginica
##   setosa         28          0         0
##   versicolor      0         21         1
##   virginica       0          3        22

31

N.0 <- 10000
N.1 <- 1000
mu.1 <- 1
mu.0 <- -1
var.1 <- 1
var.0 <- 1
x <- rnorm(N.0, mu.0, var.0)  # 病気でない人
y <- rnorm(N.1, mu.1, var.1)  # 病気の人
plot(1:1, 1:1, xlim = c(0, 1), ylim = c(0, 1),
     xlab = "False Positive", ylab = "True Positive",
     main = "ROC曲線", type = "n")
theta.seq <- exp(seq(-10, 100, 0.1))
U <- NULL
V <- NULL
for (theta in theta.seq) {
  # 病気でない人を病気とみなす
  u <- sum(dnorm(x, mu.1, var.1) / dnorm(x, mu.0, var.0) > theta) / N.0
  # 病気の人を病気とみなす
  v <- sum(dnorm(y, mu.1, var.1) / dnorm(y, mu.0, var.0) > theta) / N.1
  U <- c(U, u)
  V <- c(V, v)
}
lines(U, V)
M <- length(theta.seq) - 1
AUC <- 0
for (i in 1:M)
  AUC <- AUC + abs(U[i + 1] - U[i]) * V[i]
text(0.5, 0.5, paste("AUC = ", AUC))

略解 (Rプログラム)

第2章分類

20

22

24

26

28

29

30

31

略解 (Rプログラム)

第2章 分類

20

22

24

26

28

29

30

31

第2章分類