> head(df1)
u v
1 0.67831670 0.66337274
2 0.07438501 0.07192181
3 0.01026384 0.34820866
> head(uu)
[,1] [,2]
[1,] 0.67831670 0.66337274
[2,] 0.07438501 0.07192181
[3,] 0.01026384 0.34820866
library(ggplot2)
library(ggExtra)
library(mvtnorm)
curve(pnorm,-4,4,col=4,lty=1,lwd=4)
mtext("標準正規分布(累積分布)", col = "black", adj = 0)
text(-2.0, 0.7, expression(F(x) == paste(frac(1, sqrt(2 * pi)),
" ", integral(e^(-x^2/2) * dx, -infinity, x))), cex = 1)
arrows(1,-1,1,pnorm(1),col="red",lwd=2)
arrows(1,pnorm(1),-4.5,pnorm(1),col="green",lwd=2)
text(1.2, 0,"1.0")
text(-3.8,0.9,"0.8413")
pnorm(1)
qnorm(0.8413)
pnorm(qnorm(0.8413))
x <- seq(-4, 4, length = 401)
z<-seq(0,1,0.01)
qn<-qnorm(z)
plot(z,qn,type="l", main = "標準正規分布の逆関数 qnorm",lwd=5,col=5,
xlab= "確率点(P)",
ylab="z値")
text(0.2, 0.4, expression(F^{-1} *(z)), cex =2)
arrows(0.841,-3,0.841,1,col="red",lwd=2)
arrows(0.841,1,-0.5,1,col="green",lwd=2)
text(0, 1.2,"1.0")
text(0.75,-2,"0.8413")
set.seed(2020)
sigm<-matrix(c(1,0.7,0.7,1),byrow=T,ncol=2,nrow=2)
#平均値0 分散1 相関係数0.7の乱数2系列
samdata<-rmvnorm(n=1000,mean=c(0,0),sigma=sigm)
head(samdata)
#発生させた乱数の平均値と分散行列
colMeans(samdata)
var(samdata)
#相関係数0.7の2変数の散布図
plot(samdata[,1],samdata[,2],xlab="x",ylab="y",pch="*",
cex=1,main="同時分布(正規分布)")
#作図の都合上、matrixをdataframeに変換
dfn<-data.frame(x=samdata[,1],y=samdata[,2])
p02 <- ggplot(dfn, aes(x, y)) +geom_point() +
labs(title = "散布図)")
p02
##相関係数0.7の2変数標準正規分布の確率密度関数(PDF)
x <- y <- seq( -4, to=4, length=30)
grid <- expand.grid( x=x, y=y)
pdfgrid <- dmvnorm( grid, mean=c(0, 0), sigma=sigm)
zpdf <- array( pdfgrid, c(30, 30))
persp(x, y, zpdf, col=gray(0.9), zlab="",tick='detailed',
theta=15, phi=40,main="2次元標準正規同時密度")
##相関係数0.7の[2変数正規分布の(累積)同時分布関数(CDF)
x <- y <- seq( -4, to=4, length=30)
grid <- expand.grid( x=x, y=y)
cdfgrid <- apply(grid, 1, function(x) {
pmvnorm( upper=x, mean=c(0, 0), sigma=sigm) } )
zcdf <- array( cdfgrid, c(30, 30))
persp(x, y, zcdf, col=gray(0.9), zlab="",
tick='detailed', theta=15, phi=10,zlim=c(0,1),
main="2次元標準正規同時分布")
#確率積分変換し、同時正規分布関数に代入して、
#周辺分布が一様分布のコピュラを求める
samtocopula<-pnorm(samdata)
head(samtocopula)
plot(samtocopula,xlab="u(一様分布)",
ylab="v(一様分布)" ,
pch="*",cex=1.5,main="正規コピュラ" )
df1<-data.frame(u=samtocopula[,1],v=samtocopula[,2])
pf2 <- ggplot(df1, aes(u, v)) +geom_point()+
labs(title = "正規コピュラ")
ggMarginal(pf2, type = "histogram")
#コピュラのライブラリーcopulaを使う
library(copula)
#相関係数を0.7にセットして正規コピュラの乱数を生成し
#グラフを描く
norm.cop <- normalCopula(0.7)
set.seed(2020)
uu <- rCopula(1000, norm.cop)
plot(uu)
#生成した乱数df1とuuは等しいものになっている。
head(df1)
head(uu)
persp (norm.cop, dCopula,main="PDF",xlab="u1",ylab="u2",theta=30)
persp (norm.cop, pCopula,main="CDF",xlab="u1",ylab="u2",theta=30)
クレイトンコピュラ(Clayton copula)の自由研究