# The effects of Means and Variance on distribution
 # Written for Advanced Statistics, Econ. NTU
 # Jin-Lung Lin, March 2007
set.seed(12345678);
nobs=500;
n0=10;
p0=.25


#Binomial distribution


op0=par(mfcol=c(2,2))
curve(dbinom(x,size=n0,prob=p0),0,20,n=21,type="s")
mtext(paste("Density of Bin(",n0,",",p0,")"),side=3) 
curve(pbinom(x,size=n0,prob=p0),0,20,n=21,type="s",col="magenta")
mtext(paste("Distribution of Bin(",n0,",",p0,")"),side=3) 
curve(dnorm(x),-5,5)
mtext(paste("Density of N(0,1)"))
curve(pnorm(x),-5,5)
mtext(paste("Distribution of N(0,1)"))
par(op0)





x=c(-100,-10,-2.58,-1.96,-1.6,0,1.6,1.96,2.58,10,100)
pnorm(x)

p=c(0.005,0.025,0.05,0.5,0.95,0.975,0.995)
qnorm(p)

dnorm(x)

#   0.3989423 = 1/sqrt(2*pi) 

mu=matrix(0,3,3)
mu[1,]=-1
mu[2,]=0
mu[3,]=1

sig=matrix(0,3,3)
sig[,1]=.5
sig[,2]=1
sig[,3]=2

x0=rep(0,101)
y0=seq(0,1,by=0.01)

op=par(mfrow=c(3,3))
curve(dnorm(x,mean=mu[1,1],sd=sig[1,1]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[1,1], "sig=", sig[1,1]),side=3)

curve(dnorm(x,mean=mu[1,2],sd=sig[1,2]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[1,2], "sig=", sig[1,2]),side=3)

curve(dnorm(x,mean=mu[1,3],sd=sig[1,3]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[1,3], "sig=", sig[1,3]),side=3)

curve(dnorm(x,mean=mu[2,1],sd=sig[2,1]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[2,1], "sig=", sig[2,1]),side=3)

curve(dnorm(x,mean=mu[2,2],sd=sig[2,2]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[2,2], "sig=", sig[2,2]),side=3)

curve(dnorm(x,mean=mu[2,3],sd=sig[2,3]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[2,3], "sig=", sig[2,3]),side=3)

curve(dnorm(x,mean=mu[3,1],sd=sig[3,1]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[3,1], "sig=", sig[3,1]),side=3)

curve(dnorm(x,mean=mu[3,2],sd=sig[3,2]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[3,2], "sig=", sig[3,2]),side=3)

curve(dnorm(x,mean=mu[3,3],sd=sig[3,3]),-4,4,)
lines(x0,y0,col="magenta")
mtext(paste("mu=",mu[3,3], "sig=", sig[3,3]),side=3)
par(op)



x=rnorm(100)
y=seq(1,100)
op2=par(mfcol=c(4,2))
plot(x,y,type="p",col="red")
hist(x,prob=TRUE,ylim=c(0,1))
curve(dnorm(x),-5,5,add=T,col="magenta")
mtext(paste("Drawing from Normal(",0,1,")")   )


x=rnorm(100,0,2)
y=seq(1,100)
plot(x,y,type="p",col="red")
hist(x,prob=TRUE,ylim=c(0,1))
curve(dnorm(x,sd=2),-5,5,add=T,col="magenta")
mtext(paste("Drawing from Normal(",0,2,")")   )

x=rnorm(100,1,1)
y=seq(1,100)
plot(x,y,type="p",col="red")
hist(x,prob=TRUE,ylim=c(0,1))
curve(dnorm(x,mean=1),-5,5,add=T,col="magenta")
mtext(paste("Drawing from Normal(",1,1,")")   )



 
op3=par(mfrow=c(3,1))
curve(dnorm(x),-10,10,col=1)
mtext("Normal desnity with mean = 0, sd =1") 
curve(dnorm(x,sd=2),-10,10,col=12)
mtext("Normal desnity with mean = 0, sd =2") 
curve(dnorm(x,sd=0.5),-10,10,col=3)
mtext("Normal desnity with mean = 0, sd =0.5") 
par(op3)


curve(dnorm(x),-10,10,col=1,ylim=c(0,1))
curve(dnorm(x,sd=2),-10,10,col=2,add=T,ylim=c(0,1))
curve(dnorm(x,sd=0.5),-10,10,col=3,add=T,ylim=c(0,1))
legend(1,legend=c("N(0,1)","N(0,4)","N(0,1/4)"),lty=c(1,1,1),col=1:3)

# Gamma distributions

a=matrix(0,4,3)
b=matrix(0,4,3)
b[,1]=1
b[,2]=2
b[,3]=4
a[1,]=.5
a[2,]=1
a[3,]=1.5
a[4,]=2
op5=par(mfrow=c(4,3))
curve(dgamma(x,shape=a[1,1],scale=b[1,1]),0,15,col="magenta")
mtext(paste("Gamma(",a[1,1],",",b[1,1],")")  )
curve(dgamma(x,shape=a[1,2],scale=b[1,2]),0,15,col="magenta")
mtext(paste("Gamma(",a[1,2],",",b[1,2],")")  )
curve(dgamma(x,shape=a[1,3],scale=b[1,3]),0,15,col="magenta")
mtext(paste("Gamma(",a[1,3],",",b[1,3],")")  )
curve(dgamma(x,shape=a[2,1],scale=b[2,1]),0,15,col="magenta")
mtext(paste("Gamma(",a[2,1],",",b[2,1],")")  )
curve(dgamma(x,shape=a[2,2],scale=b[2,2]),0,15,col="magenta")
mtext(paste("Gamma(",a[2,2],",",b[2,2],")")  )
curve(dgamma(x,shape=a[2,3],scale=b[2,3]),0,15,col="magenta")
mtext(paste("Gamma(",a[2,3],",",b[2,3],")")  )
curve(dgamma(x,shape=a[3,1],scale=b[3,1]),0,15,col="magenta")
mtext(paste("Gamma(",a[3,1],",",b[3,1],")")  )
curve(dgamma(x,shape=a[3,2],scale=b[3,2]),0,15,col="magenta")
mtext(paste("Gamma(",a[3,2],",",b[3,2],")")  )
curve(dgamma(x,shape=a[3,3],scale=b[3,3]),0,15,col="magenta")
mtext(paste("Gamma(",a[3,3],",",b[3,3],")")  )
curve(dgamma(x,shape=a[4,1],scale=b[4,1]),0,15,col="magenta")
mtext(paste("Gamma(",a[4,1],",",b[4,1],")")  )
curve(dgamma(x,shape=a[4,2],scale=b[4,2]),0,15,col="magenta")
mtext(paste("Gamma(",a[4,2],",",b[4,2],")")  )
curve(dgamma(x,shape=a[4,3],scale=b[4,3]),0,15,col="magenta")
mtext(paste("Gamma(",a[4,3],",",b[4,3],")")  )


# degree of freedom in Chisquare
x0=rep(0,101)
y0=seq(0,1,by=0.01)
c=sqrt(18)
curve(dnorm(x,mean=9,sd=c),0,20,col=1)
curve(dchisq(x,df=1),0,20,col=2,add=T)
curve(dchisq(x,df=9),0,20,col=3,add=T)
lines(x0,y0)
legend("topright",legend=c("N(9,18)","Chi(1)","Chi(9)"),col=1:3,lty=c(1,1,1))

x0=rep(0,101)
y0=seq(0,1,by=0.01)
c=sqrt(40)
curve(dnorm(x,mean=20,sd=c),0,40,col=1)
curve(dchisq(x,df=1),0,40,col=2,add=T)
curve(dchisq(x,df=20),0,40,col=3,add=T)
lines(x0,y0)
legend("topright",legend=c("N(20,40)","Chi(1)","Chi(20)"),col=1:3,lty=c(1,1,1))


# Normal and Chisquare
  z=rnorm(100)
  y=z^2
  op6=par(mfrow=c(3,1))
  plot(z,y,type="h")
  hist(z,xlim=c(-10,10),prob=T)
  hist(y,xlim=c(0,10),prob=TRUE)
  curve(dchisq(x,df=1),0,10,add=T,col="magenta") 
  par(op6)

# Beta distribution
  a1=matrix(0,4,4)
  a2=matrix(0,4,4)
  a1[,1]=.5
  a1[,2]=1
  a1[,3]=2
  a1[,4]=3
  a2[1,]=3
  a2[2,]=2
  a2[3,]=1
  a2[4,]=.5
  
  op7=par(mfrow=c(4,4))
  curve(dbeta(x,shape1=a1[1,1],shape2=a2[1,1]),0,1,col="magenta")
  mtext(paste("Beta(",a1[1,1],",",a2[1,1],")")  )
  curve(dbeta(x,shape1=a1[1,2],shape2=a2[1,2]),0,1,col="magenta")
  mtext(paste("Beta(",a1[1,2],",",a2[1,2],")")  )
  curve(dbeta(x,shape1=a1[1,3],shape2=a2[1,3]),0,1,col="magenta")
  mtext(paste("Beta(",a1[1,3],",",a2[1,3],")")  )
  curve(dbeta(x,shape1=a1[1,4],shape2=a2[1,4]),0,1,col="magenta")
  mtext(paste("Beta(",a1[1,4],",",a2[1,4],")")  )
  
  curve(dbeta(x,shape1=a1[2,1],shape2=a2[2,1]),0,1,col="magenta")
  mtext(paste("Beta(",a1[2,1],",",a2[2,1],")")  )
  curve(dbeta(x,shape1=a1[1,2],shape2=a2[1,2]),0,1,col="magenta")
  mtext(paste("Beta(",a1[2,2],",",a2[2,2],")")  )
  curve(dbeta(x,shape1=a1[2,3],shape2=a2[2,3]),0,1,col="magenta")
  mtext(paste("Beta(",a1[2,3],",",a2[2,3],")")  )
  curve(dbeta(x,shape1=a1[2,4],shape2=a2[2,4]),0,1,col="magenta")
  mtext(paste("Beta(",a1[2,4],",",a2[2,4],")")  )
  
  curve(dbeta(x,shape1=a1[3,1],shape2=a2[3,1]),0,1,col="magenta")
  mtext(paste("Beta(",a1[3,1],",",a2[3,1],")")  )
  curve(dbeta(x,shape1=a1[3,2],shape2=a2[3,2]),0,1,col="magenta")
  mtext(paste("Beta(",a1[3,2],",",a2[3,2],")")  )
  curve(dbeta(x,shape1=a1[3,3],shape2=a2[3,3]),0,1,col="magenta")
  mtext(paste("Beta(",a1[3,3],",",a2[3,3],")")  )
  curve(dbeta(x,shape1=a1[3,4],shape2=a2[3,4]),0,1,col="magenta")
  mtext(paste("Beta(",a1[3,4],",",a2[3,4],")")  )
  
  curve(dbeta(x,shape1=a1[4,1],shape2=a2[4,1]),0,1,col="magenta")
  mtext(paste("Beta(",a1[4,1],",",a2[4,1],")")  )
  curve(dbeta(x,shape1=a1[4,2],shape2=a2[4,2]),0,1,col="magenta")
  mtext(paste("Beta(",a1[4,2],",",a2[4,2],")")  )
  curve(dbeta(x,shape1=a1[4,3],shape2=a2[4,3]),0,1,col="magenta")
  mtext(paste("Beta(",a1[4,3],",",a2[4,3],")")  )
  curve(dbeta(x,shape1=a1[4,4],shape2=a2[4,4]),0,1,col="magenta")
  mtext(paste("Beta(",a1[4,4],",",a2[4,4],")")  )
  par(op7)