#

# R example 3 for normal density, 2/3/06, by Yijun Zuo

#

#

#Splus build-in distribution functions include

#norm, t, f, chisq, exp, mvnorm, etc.

#

#--------------(1)-- One-dimensional-------------------------

#probability, quantile, and density of norm:

#pnorm,qnorm,dnorm (pt, qt, dt, etc.)

pnorm(0)

qnorm(0.5)

pnorm(1)

qnorm(pnorm(1))

dnorm(0)

1/(2*pi)^0.5

#------density plot--------------

yn=seq(-5,5,0.005)

plot(yn,dnorm(yn),type="l",xlab="",ylab="",main="Normal density curve",col="blue")

# 

#generating a sample from a distribution

#rnorm(100),rt(75,2),rf(100,2,3),runif(100,2,3)

#

xn=rnorm(100000)

mean(xn)

var(xn)

#

#construct density curve from the data

#

plot(density(xn),type="b",pch=".")

plot(density(xn),type="b",pch=".", main="Density curve of 1000 standard normal points",col="red")

points(yn,dnorm(yn),type="l",xlab="",ylab="",main="Normal density curve",col="blue")

#

#----------------(2)----Multi-dimensional-----------

#

#-----------load mvtnorm package first---------

#pmvnorm, qmvnorm, dmvnorm, rmvnorm

pmvnorm(c(0,0))[1]

qmvnorm(0.25,sigm=diag(2),tail="lower.tail")$quantile

dmvnorm(c(0,0))

1/(2*pi)

xn2D=rmvnorm(200,mean=c(0,0),sigma=matrix(c(1,.5,.5,1),2))

plot(xn2D,xlab="x",ylab="y",xlim=c(-4,4),ylim=c(-4,4), pch=18,col="blue")

scatterplot3d(cbind(xn2D,dmvnorm(xn2D)),color=c("red"),pch=18,highlight.3d=TRUE,

xlab="x", ylab="y", zlab="density",col.axis=c("blue"),col.grid=c("gray"))

xn2d=seq(-5,5,length=50)

persp(xn2d,xn2d,xlab="x", ylab="y", zlab="density",outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/2)/(2*pi)}))

#persp(xn2d,xn2d,outer(xn2d,xn2d,dmvnorm(c(x,y),sigma=diag(length(2)))),xlab="x",ylab="y",zlab="density")

title("Density of bivariate standard norm distribution")

contour(xn2d,xn2d,outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/2)/(2*pi)}))

contour(xn2d,xn2d,outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/2)/(2*pi)}),nlevels=3)

contour(xn2d,xn2d,outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/2)/(2*pi)}),nlevels=3)

contour(xn2d,xn2d,outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/8)/(8*pi)}),nlevels=3)

persp(xn2d,xn2d,xlab="x", ylab="y", zlab="density",outer(xn2d,xn2d,function(x,y){exp(-(x^2+y^2)/8)/(8*pi)}))

title("Density of bivariate norm distribution")