#Univariate normal distribution dnormal=function(x=0,mean=0,sd=1){ p=1/(sqrt(2*pi)*sd)*exp(-(x-mean)^2/(2*sd^2)) return(p) } x=c(95,95,95,95) mean=c(100,100,85,85) sd=c(1,2,5,10) dnormal(x,mean,sd) dnorm(x,mean,sd) dnormal(10,10,1) #Multiple variate normal distribution dmnormal=function(x=0,mean=0,V=NULL){ n=length(x) p=1/(sqrt(2*pi)^n*sqrt(det(V)))*exp(-t(x-mean)%*%solve(V)%*%(x-mean)/2) return(p) } x=matrix(c(100,95,70),3,1) mean=rep(90,3) K=matrix(c(1,.75,.25,.75,1,.25,.25,.25,1),3,3) va=95 ve=5 V=2*K*va+ve dmnormal(x,mean,V) va=5 ve=95 V=2*K*va+ve dmnormal(x,mean,V) va=50 ve=50 V=2*K*va+ve dmnormal(x,mean,V) p=c(6.897989e-06, 6.302263e-17, 3.255017e-06) output=cbind(K,x,p) par(mfrow=c(4,1)) x=rnorm(10000,100,1) plot(density(x),xlim=c(60,105)) x=rnorm(10000,100,2) plot(density(x),xlim=c(60,105)) x=rnorm(10000,85,5) plot(density(x),xlim=c(60,105)) x=rnorm(10000,85,10) plot(density(x),xlim=c(60,105))