library(mvtnorm) #多変量正規分布パッケージ
library(scatterplot3d) #三次元描画パッケージ

sigma.zero <- matrix(c(1,0,0,1), ncol=2) #分散共分散行列（無相関）
sigma.positive <- matrix(c(1,0.6,0.6,1), ncol=2) #分散共分散行列（正相関）
sigma.negative <- matrix(c(1,-0.6,-0.6,1), ncol=2) #分散共分散行列（負相関）

x100 <- rmvnorm(n=100, mean=c(0,0), sigma=sigma.zero) #乱数100個
hist(x100[ ,1]) #一次元ヒストグラム（変量１）

hist(x100[ ,2]) #一次元ヒストグラム（変量２）

plot(x100) #二次元散布図

scatterplot3d(x100[,1], x100[,2],
  dmvnorm(x100, mean=c(0,0), sigma=sigma.zero),
  highlight=TRUE) #三次元散布図

x1000 <- rmvnorm(n=1000, mean=c(0,0), sigma=sigma.zero) #乱数1000個
hist(x1000[ ,1]) #一次元ヒストグラム（変量１）

hist(x1000[ ,2]) #一次元ヒストグラム（変量２）

plot(x1000) #二次元散布図

scatterplot3d(x1000[,1], x1000[,2],
  dmvnorm(x1000, mean=c(0,0), sigma=sigma.zero),
  highlight=TRUE) #三次元散布図

x10000 <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero) #乱数10000個
hist(x10000[ ,1]) #一次元ヒストグラム（変量１）

hist(x10000[ ,2]) #一次元ヒストグラム（変量２）

plot(x10000) #二次元散布図

scatterplot3d(x10000[,1], x10000[,2],
  dmvnorm(x10000, mean=c(0,0), sigma=sigma.zero),
  highlight=TRUE) #三次元散布図

x10000posi <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.positive) #正相関
hist(x10000posi[ ,1]) #一次元ヒストグラム（変量１）

hist(x10000posi[ ,2]) #一次元ヒストグラム（変量２）

plot(x10000posi) #二次元散布図

scatterplot3d(x10000posi[,1], x10000posi[,2],
  dmvnorm(x10000posi, mean=c(0,0), sigma=sigma.positive),
  highlight=TRUE) #三次元散布図

x10000nega <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.negative) #負相関
hist(x10000nega[ ,1]) #一次元ヒストグラム（変量１）

hist(x10000nega[ ,2]) #一次元ヒストグラム（変量２）

plot(x10000nega) #二次元散布図

scatterplot3d(x10000nega[,1], x10000nega[,2],
  dmvnorm(x10000nega, mean=c(0,0), sigma=sigma.negative),
  highlight=TRUE) #三次元散布図

x1 <- seq(-3, 3, length=50)  # 変量x1の定義域 -3≦x1≦3
x2 <- seq(-3, 3, length=50)  # 変量x2の定義域 -3≦x2≦3

# 分散共分散行列 sigma.zero の密度関数
f.zero <- function(x1,x2) { 
    dmvnorm(matrix(c(x1,x2), ncol=2), 
    mean=c(0,0), sigma=sigma.zero) }
# x1とx2の定義域の外積に対する密度関数f.zeroの値域
z <- outer(x1, x2, f.zero) 
# z に関する条件
z[is.na(z)] <- 1  
#グラフィクスの環境設定
op <- par(bg = "white")  
# 密度関数の3Dグラフィクス描画
persp(x1, x2, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")  

# 分散共分散行列の固有値と固有ベクトル
eigen(sigma.zero)

## eigen() decomposition
## $values
## [1] 1 1
## 
## $vectors
##      [,1] [,2]
## [1,]    0   -1
## [2,]    1    0

# 分散共分散行列 sigma.positive の密度関数
f.posi <- function(x1,x2) { 
    dmvnorm(matrix(c(x1,x2), ncol=2), 
    mean=c(0,0), sigma=sigma.positive) }
z <- outer(x1, x2, f.posi) 
z[is.na(z)] <- 1; op <- par(bg = "white")
# 密度関数の3Dグラフィクス描画
persp(x1, x2, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")  

# 分散共分散行列の固有値と固有ベクトル
eigen(sigma.positive)

## eigen() decomposition
## $values
## [1] 1.6 0.4
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

# 分散共分散行列 sigma.negative の密度関数
f.nega <- function(x1,x2) { 
    dmvnorm(matrix(c(x1,x2), ncol=2), 
    mean=c(0,0), sigma=sigma.negative) }
z <- outer(x1, x2, f.nega) 
z[is.na(z)] <- 1; op <- par(bg = "white")
# 密度関数の3Dグラフィクス描画
persp(x1, x2, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")  

# 分散共分散行列の固有値と固有ベクトル
eigen(sigma.negative)

## eigen() decomposition
## $values
## [1] 1.6 0.4
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068