########################################################
###### Numerical Methods & Simulation ##################
########################################################

### class 1

## creat a function 

#

x<--2:2

x<-seq(-2,2,by=0.1)
x
y<-x^2
plot(x,y,type="l")


## curve

curve(x^2,-10,10)
curve(x^3,-10,10)
curve(y^3,-10,10)  ## error

#############################################
# class 2
##############################################
## problem 1: plot a function in R by "curve"
# 1: f(x)= x*log(x)     x>0
#2: f(x)=x*sin(2*pi*x)    pi is 3.14,
#3: f(x)= log(x)/(1+x)   x>0
#4: f(x)=1/sqrt(2*pi)*exp(-x^2/2)      f(x) is density of the standard normal distribution


# 3: f(x)= log(x)/(1+x)   x>0

curve(log10(x)/(1+x),0,10)
abline(v=c(2,4),lty=2,col=2)

# create a function

f<-function(x) {
	y<-log10(x)/(1+x)
	return(y)
}

f(2)
f(4)
f(10)
f(50)

curve(f(x),0,10)


f1<-function(x,theta=.5) {
	y<-log10(x)/(1+theta*x)
	return(y)
}

f1(x=2,theta=5)
f1(2,5)
f1(5,10)
f1(5)

curve(f1(x,theta=1),0,10)
curve(f1(x,theta=2),0,10,add=TRUE,lty=2,col=2)
curve(f1(x,theta=5),0,10,add=TRUE,lty=3,col=3)


## normal density plots with different mean

N1<-function(x,mu=0,sigma=1){
	1/sqrt(2*pi*sigma^2)*exp(-(x-mu)^2/(2*sigma^2))
}

N1(2,mu=2,sigma=2)
curve(N1(x,0,1),-5,5)
curve(N1(x,1,1),-5,5,add=T,col=2,lty=2)

#############################################
# class 3
##############################################
########## screen Bedore
N1<-function(x,mu,sigma){
  1/sqrt(2*pi*sigma^2)*exp(-(x-mu)^2/(2*sigma^2))
}
N1(2,2,2)
curve(N1(x,0,1),-5,5,ylab = "f(x)",ylim=c(0,0.8))
curve(N1(x,0,2),-5,5,add=T,col=2,lty=2)
curve(N1(x,0,4),-5,5,add=T,col=3,lty=3)
curve(N1(x,0,0.5),-5,5,add=T,col=4,lty=4)
legend("topleft",legend = c("N(0,1)","N(0,2)","N(0,4)","N(0,0.5)"),col=1:4,lty=1:4)
?legend
##Fisher##
?distribution
##############################
F1<-function(x,n1,n2){
 gamma((n1+n2)/2)/(gamma(n1/2)*gamma(n2/2))*((n1/n2)^(n1/2))*x^((n1/2)-1)/(1+n1/n2*x)^((n1+n2)/2)
   }
F1(1,1,1)
curve(F1(x,1,1),0,10)   
curve(F1(x,2,3),add = T,col=2,lty=2)
curve(F1(x,4,5),add = T,col=3,lty=3)
###########################################
df(1,1,1)
curve(df(x,1,1),0,10)   
curve(df(x,1,3),add = T,col=2,lty=2)
curve(df(x,1,5),add = T,col=3,lty=3)
curve(df(x,1,15),add = T,col=4,lty=4)

##Chi Squared##
#######################################
### Error
A1<-function(x,k){
  x^(k/2)*exp(-(x/2))/(2^(k/2)*gamma(k/2))
}
A1(1,1)
curve(A1(x,1))
curve(A1(x,2),add = T,col=2,lty=2)
#####################################

dchisq(1,1)
##curve(dchisq(x,1),0,10)
curve(dchisq(x,2),0,30,col=2,lty=2)
curve(dchisq(x,3),0,30,add = T,col=3,lty=3)
curve(dchisq(x,5),0,30,add = T,col=4,lty=4)
curve(dchisq(x,10),0,30,add = T,col=5,lty=5)
curve(dchisq(x,15),0,30,add = T,col=6,lty=6)
?distribution
 
#############################################
# class 4
##############################################
# create likelihood functions by n random samples

#Example1: f(x)=(1+a*x)/2     -1<x<1    -1<a<1


### for? loops in R ?

#?for
for (i in vector) {
   ...
}


for(i in 1:10) {
	print(i)
}

w=0
for(j in c(2,5,12,-5)){
	print(j)
	w<-j+w
}
w

z<-c(3,5,7,13,-4)

z[1]
z[2]
z[3]
z[4]
z[5]
z[6]

m<-matrix(c(1,2,3,4,5,6),nrow=3,ncol=2)
m

m[1,1]
m[2,2]

m[1,]

m[,2]

x=0
y=0
for(k in 1:10){
	x[k]<-k^3
	y[k]<-k*(k+1)
}

x
y

# f(x)=e^(-k/2)  k= 1,2,3,...,100
f=0
for(k in 1:100){
	f[k]<-exp(-k/2)
}
f

#example: for

# f(x)=x*(x+1)*r  x= 1,2,3,...,10   0<r<1

# plot function prod(f(x)) by r?

# home work



# g(x)=exp(-x/b)  x>0 , b>0 
x1<-c(0.1,0.2,0.5,0.8,0.4,0.5,0.6,0.3,0.6,0.8) # 10 random sample from g(x)

# plot g(x) by b?

b<-0.5
g<-0
for(i in 1:10){
	g[i]<-exp(-x[i]/b)
}

g

# create function g by b

g<-function(b,x){
	g<-0
    for(i in 1:10){
	 g[i]<-exp(-x[i]/b)
    }
 return(g)
}

g(.5,x1)

g(0.6,x1)


f<-function(b,x1){
	f<-matrix(0,ncol=length(x1),nrow=length(b))
for(j in 1:length(b)){
	f[j,]<-g(b[j],x1)
}
f<-apply(f,1,prod)
return(f)
}

f(1:100,x1)
curve(f(x,x1),0,200)
curve(f(x,x1),0,1)


###################################################
# simple function by log
# logarithm
# edit function by logarithm
# prod convert to sum
# log(g(x))=-x/b  then logf(x1,...,x10)=sum(-xi/b)

 logf<-function(b,x1){
	f<-0
for(j in 1:length(b)){
	f[j]<-sum(-x1/b[j])
}
return(f)
}

b=c(0.5,0.6)

logf(b,x1)
f(b,x1)

par(mfrow=c(1,2))
curve(logf(x,x1),0,0.1)
curve(f(x,x1),0,0.1)

##################
### problem 
# generating 10 random samples from normal(0,1)
# likelihood function , mu=?: f(mu,x2)
n=10
set.seed(134)
x2<-rnorm(n)
x2

# create likelihood function by mu and draw it
n<-10

f<-function(mu){
  set.seed(134)
  x<-rnorm(n)
  f<-matrix(0,ncol=length(x),nrow = length(mu))
  for (i in 1:length(mu)) {
    f[i,]<-dnorm(x,mu[i])
  }
  f<-apply(f,1,prod)
  return(f)
}
f(c(0,1,2))
curve(f(x),-2,2)

curve(dnorm(x,2),-3,6)

set.seed(134)
x<-rnorm(n)
mean(x)
abline(v = c(mean(x),0),lty = 2,col = 2)
#############################################
# class 5
##############################################
# likelihood function of Poisson distribution
#f(x,lambda)=exp(-lambda)*lambda^(x)/factorial(x)
set.seed(143)
y<-rpois(10,4)
y


Pois<-function(lambda,y){
	f<-0
	for(i in 1:length(lambda)){
		f[i]<-prod(exp(-lambda[i])*lambda[i]^(y)/factorial(y))
	}
	return(f)
}
lambda<-c(1,2,3)
Pois(lambda,y)

curve(Pois(x,y),1,10)

### log poisson


log.Pois<-function(lambda,y){
	f<-0
	for(i in 1:length(lambda)){
		f[i]<-sum(-lambda[i]+y*log(lambda[i])-log(factorial(y)))
	}
	return(f)
}
lambda<-c(1,2,3)
log.Pois(lambda,y)

par(mfrow=c(1,2))
curve(Pois(x,y),1,10)
curve(log.Pois(x,y),1,10)





# minimize f
optimize(function(b) f(b,x1), interval=c(0,100))
optim(.5,function(b) f(b,x1),method="Brent",lower=0,upper=100)
optim(.5,function(b) f(b,x1),method="BFGS")
optim(.5,function(b) f(b,x1),method="Nelder-Mead")
optim(.5,function(b) f(b,x1),method="CG")
optim(.5,function(b) f(b,x1),method="L-BFGS-B",lower=0,upper=100) 
optim(.5,function(b) f(b,x1),method="SANN")

# g(x)=exp(-x/b)  x>0 , b>0 
x1<-c(0.1,0.2,0.5,0.8,0.4,0.5,0.6,0.3,0.6,0.8) # 10 random sample from g(x)


# simple function by log
# logarithm
# edit function by logarithm
# prod convert to sum
# log(g(x))=-x/b  then logf(x1,...,x10)=sum(-xi/b)

 logf<-function(b,x1){
	f<-0
for(j in 1:length(b)){
	f[j]<-sum(-x1/b[j])
}
return(f)
}
b=c(0.5,.6)
logf(b,x1)

curve(logf(x,x1),0,1)
# minimize logf
optimize(function(b) logf(b,x1), interval=c(0,10))

optim(.5,function(b) logf(b,x1),method="Brent",lower=0,upper=1)

optim(.5,function(b) logf(b,x1),method="BFGS")

optim(.5,function(b) logf(b,x1),method="Nelder-Mead")

optim(.5,function(b) logf(b,x1),method="CG")

optim(.5,function(b) logf(b,x1),method="L-BFGS-B",lower=0,upper=100) 
optim(.5,function(b) logf(b,x1),method="SANN")

Min<-optimize(function(b) logf(b,x1), interval=c(0,10))
Min$minimum
curve(logf(x,x1),0,.01)
abline(v=Min$minimum,lty=2,col=2)

#############################################
# class 6
##############################################

### problem 3 distributions ###

### code khodavirdi:
###likelihood function of exponential distribution
set.seed(143)
y<-rexp(10,5)
y
e<-function(lambda,x){
  f<-0
  for(i in 1:length(lambda)){
    f[i]<-prod(lambda[i]*exp(-lambda[i]*x))
  }
  return(f)
}
lambda<-c(1,2)
e(lambda,y)
curve(e(x,y),0,7,col=2)

###likelihood function of chi-square distribution

set.seed(123)
y<-rchisq(100,10)
y

CHisq<-function(df,x){
  f<-0
  for(i in 1:length(df)){
    f[i]<-prod(dchisq(x,df[i]))
  }
  return(f)
}
df<-c(1,2,5)
CHisq(df,y) 
curve(CHisq(x,y),1,20)  

###likelihood function of student-t distribution
set.seed(321)
y<-rt(100,30)
y

Tstu<-function(df,x){
  f<-0
  for(i in 1:length(df)){
    f[i]<-prod(dt(x,df[i]))
}
  return(f)
}
df<-c(1,2,5)
Tstu(df,y)  
curve(Tstu(x,y),1,40) 



#############################################
# class 7
##############################################

##################
###### problem 
# create likelihood functions by n random samples
# x is c(0.41,0.91,-0.61,0.38,0.37,0.36,0.01,-0.28,-0.33,0.99,-0.35,0.1,0.31,0.75,-0.34) 

#Example1: f(x)=(1+a*x)/2     -1<x<1    -1<a<1


   
x<- c(.41,.91,-.61,.38,.37,.36,.01,-.28,-.33,.99,-.35,.1,.31,.75,-.34)

 f<- function(a){
   z<-0
   for(i in 1:length(a)){
     z[i]<-prod((1+a[i]*x)/2)
     
   }
   return(z)
 }
 

 f(c(-.5,.5,.7))
 curve(f(x),.5,1)
  abline(v = ".8",lty = 2,col = 2)
  
#### bernouli likelihood
set.seed(143)
n = 30
y<-rbinom(n,1,.5)
y
Bernuli<-function(p,x){
  f<-0
  for(i in 1:length(p)){
    f[i]<-sum(x*log(p[i])+(1-x)*log(1-p[i]))
    
  }
  return(f)
}


p<-(c(.3,.8,.1))
Bernuli(p,y)
curve(Bernuli(x,y),0,1)
abline(v = .5,lty = 2,col = "red")


#############################################
# class 8
##############################################

################################
# gamma distribution Likelihood
###############################
# f(x)= 1/(beta^alpha*Gamma(alpha))* x^(alpha-1)*e^-(x/beta)  ---  x~ G(alpha , beta)
# E(X)=alpha*beta

# beta is fixed in one (beta=1)
L.Gamma <- function(alpha){
	beta=1
	n=100
	set.seed(143)
	y<-rgamma(n,10,1)
	f<-0
	for(i in 1:length(alpha)){
	f[i]<-prod(1/(beta^alpha[i]*gamma(alpha[i]))*y^(alpha[i]-1)*exp(-y/beta))
	}
	return(f)
}
	
alpha<-c(1,4,7)
L.Gamma(alpha)

curve(L.Gamma(x),1,20)
	

# log likelihood of gamma distribution
# logf(x)=-alpha*log(beta)-log(gamma(alpha))+(alpha-1)*log(x)-x/beta

# beta is fixed in one (beta=1)
Log.Gamma <- function(alpha){
	beta=1
	n=100
	set.seed(143)
	y<-rgamma(n,10,1) # real alpha is 10
	f<-0
	for(i in 1:length(alpha)){
	f[i]<-sum(-alpha[i]*log(beta)-log(gamma(alpha[i]))+(alpha[i]-1)*log(y)-y/beta)
	}
	return(f)
}
	
alpha<-c(1,4,7)
Log.Gamma(alpha)

curve(Log.Gamma(x),1,20)

par(mfrow=c(1,2))
curve(Log.Gamma(x),5,15,xlab="alpha")
title("Log likelihood of gamma distribution")
text(9,-110,"real alpha",col=2,cex=.8)
abline(v=10,lty=2,col=2,lwd=1.5)
curve(L.Gamma(x),5,15,xlab="alpha")
title("Likelihood of gamma distribution")
text(9,5e-34,"real alpha",col=2,cex=.8)
abline(v=10,lty=2,col=2,lwd=1.5)


# minimize logf
a1<-optimize(function(alpha) (-1)*L.Gamma(alpha), interval=c(9,12))
a1
a2<-optimize(function(alpha) (-1)*Log.Gamma(alpha), interval=c(9,12))

alphahat1<-a1$minimum
alphahat2<-a2$minimum

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="Brent",lower=9,upper=11)

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="BFGS")

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="Nelder-Mead")

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="CG")

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="L-BFGS-B",lower=9,upper=11) 

optim(9.5,function(alpha) (-1)*Log.Gamma(alpha),method="SANN")


par(mfrow=c(1,2))
curve(Log.Gamma(x),5,15,xlab="alpha")
title("Log likelihood of gamma distribution")
text(9,-110,"real alpha",col=2,cex=.8)
text(11.5,-110,"alphahat",col=4,cex=.8)
abline(v=10,lty=2,col=2,lwd=1.5)
abline(v=alphahat2,,lty=3,col=4,lwd=1.5)
curve(L.Gamma(x),5,15,xlab="alpha")
title("Likelihood of gamma distribution")
text(9,5e-34,"real alpha",col=2,cex=.8)
text(11.5,5e-34,"alphahat",col=4,cex=.8)
abline(v=10,lty=2,col=2,lwd=1.5)
abline(v=alphahat1,,lty=3,col=4,lwd=1.5)

#############################################
# class 9
##############################################
## two parameters , computing likelihood

# Example: X~gamma(alpha , beta) , we have two parameters as (alpha, beta)
# log likelihood of gamma distribution
# logf(x)=-alpha*log(beta)-log(gamma(alpha))+(alpha-1)*log(x)-x/beta

Log.Gamma2<- function(alpha,beta){
	n=100
	set.seed(143)
	y<-rgamma(n,10,1) 
	f<-matrix(0,ncol=length(beta),nrow=length(alpha))
	for(i in 1:length(alpha)){
		for(j in 1:length(beta)){
		 f[i,j]<-sum(-alpha[i]*log(beta[j])-log(gamma(alpha[i]))+(alpha[i]-1)*log(y)-y/beta[j])
		}
	}
	return(f)
}
	
alpha<-c(1,4,7)
beta<-c(1,5,8,10)
Log.Gamma2(alpha,beta)

curve(Log.Gamma2(x,beta=1),1,20)

par(mfrow=c(2,2))
curve(Log.Gamma2(x,beta=1),1,20)
curve(Log.Gamma2(x,beta=2),1,20)
curve(Log.Gamma2(x,beta=4),1,20)
curve(Log.Gamma2(x,beta=.5),1,20)


#ploting surface
# by persp and image
alpha<-seq(1,25,l=20)
beta<-seq(.1,5,l=20)

f<-Log.Gamma2(alpha,beta)

persp(alpha,beta,f)

par(mfrow=c(1,2))
persp(alpha, beta, f, theta = 45, phi = 25,r=2,col = "springgreen", shade = 0.5, nticks=5, ticktype="detailed")

image(alpha, beta, f)
contour(alpha, beta, f,add=T)

### optim is general optimization
## nlm is Non-Linear Minimization by a Newton-type algorithm 
## mle library(stats4) is Maximum Likelihood Estimation
### mle is similar to optim

optim(c(1,1),function(x) (-1)*Log.Gamma2(alpha=x[1],beta=x[2]),method="Nelder-Mead")

nlm(function(x) (-1)*Log.Gamma2(alpha=x[1],beta=x[2]),c(1,1))

mle ?


#problem= beta distribution, F distribution , normal distribution

############## Beta distribution ###############
## two parameters 
Log.Beta2<- function(alpha,Beta){
  n=1000
  set.seed(143)
  y<-rbeta(n,10,1) 
  f<-matrix(0,ncol=length(Beta),nrow=length(alpha))
  for(i in 1:length(alpha)){
    for(j in 1:length(Beta)){
      f[i,j]<-sum((alpha[i]-1)*log(y)+(Beta[j]-1)*log(1-y)-log(beta(alpha[i],Beta[j])))
    }
  }
  return(f)
}

alpha<-c(1,3,5,7)
Beta<-c(2,4,6,8,10)
Log.Beta2(alpha,Beta)

curve(Log.Beta2(x,Beta=1),1,20)

par(mfrow=c(2,2))
curve(Log.Beta2(x,Beta=1),1,20)
curve(Log.Beta2(x,Beta=2),1,20)
curve(Log.Beta2(x,Beta=4),1,20)
curve(Log.Beta2(x,Beta=.5),1,20)


#ploting surface
# by persp and image
alpha<-seq(1,25,l=20)
Beta<-seq(.1,5,l=20)

f<-Log.Beta2(alpha,Beta)

persp(alpha,Beta,f)

par(mfrow=c(1,2))
persp(alpha, Beta, f, theta = 30, phi = 25,r=2,col = "springgreen", shade = 0.5, nticks=5, ticktype="detailed")

image(alpha, Beta, f)
contour(alpha, Beta, f,add=T)


optim(c(1,1),function(x) (-1)*Log.Beta2(alpha=x[1],Beta=x[2]),method="Nelder-Mead")

nlm(function(x) (-1)*Log.Beta2(alpha=x[1],Beta=x[2]),c(1,1))



#############################################
# class 10
##############################################
## function g(x) without loop (for)  by matrix function "outer"
# example: outer(x,y,FUN=function(x,y) x+y)

g_new<-function(b,x){
	f<-outer(b,x,FUN=function(b,x) exp(-x/b))
	f<-apply(f,1,prod)
	return(f)
}
b<-c(0.5,0.6,0.01,0.9)
y<-c(0.1,0.2,0.5,0.8,0.4,0.5,0.6,0.3,0.6,0.8) # 10 random sample from g(x)
g_new(b,y)

g_new(1:10,y)

curve(g_new(x,y),0,2)


## gamma(alpha,beta) : Create the Likelihood without for!!

Log.g<-function(alpha,x){
	f<-outer(alpha,x,FUN=function(alpha,x) dgamma(x,alpha,1,log=TRUE))
	f<-apply(f,1,sum)
	return(f)
}
alpha<-c(1,3,5)
set.seed(143)
n=1000
y<-rgamma(n,2,1)
Log.g(alpha,y)

curve(Log.g(x,y),0,5)
abline(v=2,lty=2,col=2)

# one line program!!!

apply(outer(c(1,3,5),rgamma(1000,2,1),FUN=function(alpha,x) dgamma(x,alpha,1,log=TRUE)),1,sum)


curve(apply(outer(x,rgamma(1000,2,1),FUN=function(alpha,y) dgamma(y,alpha,1,log=TRUE)),1,sum),0,5)

#############################################
# class 12
##############################################
######## Beta liklihood one parameters #################
Beta.lik<-function(B,x){
  f<-outer(B,x,FUN=function(B,x) dbeta(x,1,B))
  f<-apply(f,1,prod)
  return(f)
  
}

B<-c(1,3,5)
set.seed(143)
n=30
y<-rbeta(n,1,2)
Beta.lik1<-function(Beta,x){
  k<-0
  for(i in 1:length(Beta)){
    k[i]<-prod(dbeta(x,1,Beta[i]))
  }
  return(k)
}


Beta.lik(B,y)
Beta.lik1(B,y)
curve(Beta.lik(x,y),0,6)
curve(Beta.lik1(x,y),0,6,add = TRUE,col = 4,lty = 2)
abline(v=c(2,2.3),lty=2,col=c(2,3))


###### Beta log-liklihoood #######################
Log.B<-function(B,x){
  f<-outer(B,x,FUN=function(B,x) dbeta(x,B,1,log=TRUE))
  f<-apply(f,1,sum)
  return(f)
  
}

B<-c(1,3,5)
set.seed(143)
n=500
y<-rbeta(n,2,1)
Log.B(B,y)
result<-optim(1,function(x) (-1)*Log.B(B=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Log.B(x,y),0,5)
abline(v=c(2,result$par),lty=2,col=2:3)


### plot different samples
par(mfrow = c(2,2))
par(mar = rep(2,4))
set.seed(143)
n=10
y<-rbeta(n,2,1)
result<-optim(1,function(x) (-1)*Log.B(B=x,y),method="Brent",lower = .5,upper = 5)
curve(Log.B(x,y),0,5)
title("n=10")
abline(v=c(2,result$par),lty=2,col=2:3)

set.seed(143)
n=30
y<-rbeta(n,2,1)
result<-optim(1,function(x) (-1)*Log.B(B=x,y),method="Brent",lower = .5,upper = 5)
curve(Log.B(x,y),0,5)
title("n=30")
abline(v=c(2,result$par),lty=2,col=2:3)

set.seed(143)
n=100
y<-rbeta(n,2,1)
result<-optim(1,function(x) (-1)*Log.B(B=x,y),method="Brent",lower = .5,upper = 5)
curve(Log.B(x,y),0,5)
title("n=100")
abline(v=c(2,result$par),lty=2,col=2:3)


set.seed(143)
n=500
y<-rbeta(n,2,1)
result<-optim(1,function(x) (-1)*Log.B(B=x,y),method="Brent",lower = .5,upper = 5)
curve(Log.B(x,y),0,5)
title("n=500")
abline(v=c(2,result$par),lty=2,col=2:3)

#############################################
# class 13
##############################################
################### likelihood fisher  #############################
set.seed(143)
n=10
y<-rf(10,0.3,0.4)
df2<-2
fisher<-function(df1,y){
  f<-0
  for(i in 1:length(df1)){
    f[i]<-prod(df(y,df1[i],df2))
  }
  return(f)
}
df1<-c(1,2)
fisher(df1,y)
curve(fisher(x,y),0,5,col=4)

##################### likelihood fisher without for  ########################
lik.fisher1<-function(df1,x){
  f<-outer(df1,x,FUN=function(df1,x) df(x,df1,2))
  f<-apply(f,1,prod)
  return(f)
  
}
par(mfrow = c(2,2))
df1<-c(1,3,5)
set.seed(143)
n=10
y<-rf(n,1,2)
lik.fisher1(df1,y)
result<-optim(1,function(x) (-1)*lik.fisher1(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(lik.fisher1(x,y),0.1,4)
abline(v=c(1,result$par),lty=5,col=4:5)
title("n=10")

set.seed(143)
n=30
y<-rf(n,1,2)
lik.fisher1(df1,y)
result<-optim(1,function(x) (-1)*lik.fisher1(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(lik.fisher1(x,y),0.1,4)
abline(v=c(1,result$par),lty=2,col=2:3)
title("n=30")

set.seed(143)
n=100
y<-rf(n,1,2)
lik.fisher1(df1,y)
result<-optim(1,function(x) (-1)*lik.fisher1(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(lik.fisher1(x,y),0.1,4)
abline(v=c(1,result$par),lty=6,col=2:3)
title("n=100")

set.seed(143)
n=500
y<-rf(n,1,2)
lik.fisher1(df1,y)
result<-optim(1,function(x) (-1)*lik.fisher1(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(lik.fisher1(x,y),0.1,4)
abline(v=c(1,result$par),lty=4,col=2:3)
title("n=500")

######################  log likelihoode  fisher #######################
loglik.fisher<-function(df1,y){
  f<-0
  for(i in 1:length(df1)){
    f[i]<-sum(df(y,df1[i],2,log =TRUE))
  }
  return(f)
}
set.seed(143)
n=10
y<-rf(n,1,2)
df1<-c(1,2)
loglik.fisher(df1,y)
result<-optim(1,function(x) (-1)*loglik.fisher(df1=x,y),method="Brent",lower = .5,upper = 5)
result
curve(loglik.fisher(x,y),0.1,4,col=4)
abline(v=c(1,result$par),lty=4,col=2:3)

######################  log likelihood  fisher without for #######################
Loglik.fisherl<-function(df1,x){
  f<-outer(df1,x,FUN=function(df1,x) df(x,df1,2,log=TRUE))
  f<-apply(f,1,sum)
  return(f)
  
}
par(mfrow = c(2,2))
df1<-c(1,3,5)
set.seed(143)
n=10
y<-rf(n,1,2)
Loglik.fisherl(df1,y)
result<-optim(1,function(x) (-1)*Loglik.fisherl(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Loglik.fisherl(x,y),0.1,4)
abline(v=c(1,result$par),lty=4,col=2:3)
title("n=10")

set.seed(143)
n=500
y<-rf(n,1,2)
Loglik.fisherl(df1,y)
result<-optim(1,function(x) (-1)*Loglik.fisherl(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Loglik.fisherl(x,y),0.1,4)
abline(v=c(1,result$par),lty=4,col=2:3)
title("n=500")


set.seed(143)
n=1000
y<-rf(n,1,2)
Loglik.fisherl(df1,y)
result<-optim(1,function(x) (-1)*Loglik.fisherl(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Loglik.fisherl(x,y),0.1,4)
abline(v=c(1,result$par),lty=4,col=2:3)
title("n=1000")


set.seed(143)
n=5000
y<-rf(n,1,2)
Loglik.fisherl(df1,y)
result<-optim(1,function(x) (-1)*Loglik.fisherl(df1=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Loglik.fisherl(x,y),0.1,4)
abline(v=c(1,result$par),lty=4,col=2:3)
title("n=5000")

#########################


## look at the Book! (Meshkani)





#############################################
# class 13
##############################################
####norm log-liklihoood
Loglik.norm<-function(sigma,x){
  f<-outer(sigma,x,FUN=function(sigma,x) dnorm(x,1,sigma,log=TRUE))
  f<-apply(f,1,sum)
  return(f)
  
}

sigma<-c(1,3,5)
set.seed(143)
n=10
y<-rnorm(n,1,2)
Loglik.norm(sigma,y)
result<-optim(1,function(x) (-1)*Loglik.norm(sigma=x,y),method="Brent",lower = .5,upper = 5)
result
result$par
curve(Loglik.norm(x,y),0,4.5)
abline(v=c(2,result$par),lty=2,col=2:3)

par(mfrow = c(2,2))
par(mar = rep(2,4))

set.seed(143)
n=500
y<-rnorm(n,1,2)
result<-optim(1,function(x) (-1)*Loglik.norm(sigma=x,y),method="Brent",lower = .5,upper = 6)
curve(Loglik.norm(x,y),0,5)
title("n=500")
abline(v=c(2,result$par),lty=2,col=2:5)

set.seed(143)
n=1000
y<-rnorm(n,1,2)
result<-optim(1,function(x) (-1)*Loglik.norm(sigma=x,y),method="Brent",lower = .5,upper = 6)
curve(Loglik.norm(x,y),0,5)
title("n=1000")
abline(v=c(2,result$par),lty=2,col=2:5)

set.seed(143)
n=5000
y<-rnorm(n,1,2)
result<-optim(1,function(x) (-1)*Loglik.norm(sigma=x,y),method="Brent",lower = .5,upper = 6)
curve(Loglik.norm(x,y),0,5)
title("n=5000")
abline(v=c(2,result$par),lty=2,col=2:5)

set.seed(143)
n=10000
y<-rnorm(n,1,2)
result<-optim(1,function(x) (-1)*Loglik.norm(sigma=x,y),method="Brent",lower = .5,upper = 6)
curve(Loglik.norm(x,y),0,5)
title("n=10000")
abline(v=c(2,result$par),lty=2,col=2:5)

############### beta distribution two parameters

Loglik.Beta<-function(alpha,Beta,x){
  g<-0
  for(i in 1:length(alpha)){
    f<-outer(Beta,x,FUN=function(Beta,x) dbeta(x,alpha[i],Beta,log=TRUE))
    f<-apply(f,1,sum)
    g<-c(g,f)
  }
  g<-g[-1]
  return(g)
}
alpha<-1:3
Beta<-4:5
set.seed(143)
n=100
y<-rbeta(n,2,5)
Loglik.Beta(alpha,Beta,y)
result<-optim(c(1,1),function(x) (-1)*Loglik.Beta(alpha=x[1],Beta=x[2],y),method="BFGS")
result

#############################################
# class 14
##############################################
# Poisson MLE and confidence MLE and plot of MLE
# read the data 


y <- c(rep(0, 14), rep(1, 30), rep(2, 36), rep(3, 68), rep(4, 43), rep(5, 
	43), rep(6, 30), rep(7, 14), rep(8, 10), rep(9, 6), rep(10, 4), 
	rep(11, 1), rep(12, 1))
	
y<-rep(0:12,c(14,30,36,68,43,43,30,14,10,6,4,1,1))

table(y)

t(t(table(y)))

hist(y)
barplot(table(y),col="orange")
abline(v=mean(y),col="blue",lty=2,lwd=3)
mean(y)

qqplot(y,rpois(10000,mean(y)))
abline(c(0,1),col=2,lty=2)


qqplot(y,rpois(10000,5),ylim=c(0,12))
abline(c(0,1),col=2,lty=2)

?ks.test
ks.test(y,rpois(1000,mean(y)))

ks.test(y,rpois(1000,5))


n <- length(y)

Nloglike <- function(lambda,y) {
	f<- outer(y,lambda,FUN= function(y,lambda) y*log(lambda)-lambda-log(factorial(y)))
	f<-(-1)*apply(f,2,sum)
	return(f)
}
Nloglike(1:10,y)



# using MLE to find the estimated value 
result<-optim(0.5,Nloglike, y=y , method = "Brent", lower = 0, upper = 10, hessian = TRUE)
result

lambda.hat<-result$par
sd.mle<-sqrt(result$hessian^(-1))
sd.mle


lambda.hat
# install.packages("knitr")
library(knitr)
print(kable(cbind(direct=c(mean=mean(y), sd=sqrt(mean(y))), optim =c(result$par,sd.mle)),digits = 3))



# MLE function for Poission Distribution 
# and Lambda0 is starting value in optim function
MLE.Pois<-function(lambda0,y){  
	Nloglike <- function(lambda,y) {
	f<- outer(y,lambda,FUN= function(y,lambda) y*log(lambda)-lambda-log(factorial(y)))
	f<-(-1)*apply(f,2,sum)
	return(f)
}
result<-optim(lambda0,Nloglike,y=y,method = "Brent", lower = 0, upper = 10, hessian = TRUE)
return(result)
	}

Re<-MLE.Pois(0.5,y)

lambda.hat<-Re$par

# bootstraping for find distribution of MLE

B=2000   ## iteration Bootstrap
n<-length(y)
Lhat.Bootstrap<-numeric(B)
for(i in 1:B){
	ystar<-rpois(n,lambda.hat)
	Lhat.Bootstrap[i]<-MLE.Pois(0.5,ystar)$par
}


# histogram
hist(Lhat.Bootstrap, main = "sampling distribution of Bootrapped MLE for Poisson", xlab = "estimated value", breaks = 30, prob = TRUE,col="gray95")
abline(v=mean(Lhat.Bootstrap),col=2,lty=2,lwd=2)
lines(density(Lhat.Bootstrap),col="blue",lty=3)
mu<-mean(Lhat.Bootstrap)
sigma<-sd(Lhat.Bootstrap)
curve(dnorm(x,mu,sigma),col=colors()[91],add=TRUE)

### Another Version of Bootstrap 
B=2000   ## iteration Bootstrap
n<-length(y)
Lhat.Bootstrap1<-numeric(B)
for(i in 1:B){
	ystar1<-sample(y,n,replace=TRUE)
	Lhat.Bootstrap1[i]<-MLE.Pois(0.5,ystar1)$par
}


# histogram
hist(Lhat.Bootstrap1, main = "sampling distribution of Bootrapped MLE for Poisson", xlab = "estimated value", breaks = 30, prob = TRUE,col="gray95")
abline(v=mean(Lhat.Bootstrap1),col=2,lty=2,lwd=2)
lines(density(Lhat.Bootstrap1),col="blue",lty=3)
mu<-mean(Lhat.Bootstrap1)
sigma<-sd(Lhat.Bootstrap1)
curve(dnorm(x,mu,sigma),col=colors()[91],add=TRUE)

### comparing methods 
plot(density(Lhat.Bootstrap),col="blue",lty=2)
lines(density(Lhat.Bootstrap1),col="red",lty=3)
abline(v=mean(Lhat.Bootstrap),col= "blue",lty=2)
abline(v=mean(Lhat.Bootstrap1),col= "red",lty=3)


par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(Lhat.Bootstrap,col="gray70")
abline(h=mean(Lhat.Bootstrap),col="red",lty=2)

ts.plot(Lhat.Bootstrap1,col="gray70")
abline(h=mean(Lhat.Bootstrap1),col="red",lty=2)


# confidence interval by bootstraping methods

# Hessian methods
Re<-MLE.Pois(0.5,y)
Re
lambda.hat<-$Re
SE.lambda<-sqrt((Re$hessian)^(-1))
alpha=0.05
L<-lambda.hat-qnorm(1-alpha/2)*SE.lambda
L
U<-lambda.hat+qnorm(1-alpha/2)*SE.lambda
U
c(L,U)

# Bootstrap methods

conf.boot<-quantile(Lhat.Bootstrap, probs=c(0.025,0.975), names=FALSE) 
conf.boot
conf.boot1<-quantile(Lhat.Bootstrap1, probs=c(0.025,0.975), names=FALSE) 
conf.boot1

## plot confidence intervals

par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(Lhat.Bootstrap,col="gray70")
abline(h=mean(Lhat.Bootstrap),col="red",lty=2)
abline(h=c(L,U),lty=3,col="orange")



ts.plot(Lhat.Bootstrap1,col="gray70")
abline(h=mean(Lhat.Bootstrap1),col="red",lty=2)
abline(h=conf.boot,lty=3,col="orange")
abline(h=conf.boot1,lty=3,col="green")


polygon(c(1:B, rev(1:B)), c(rep(conf.boot[1],B), rev(rep(conf.boot[2],B))), col = rgb(.9,.2,.8,alpha=.2), border = NA)
##
#############################################
# class 17
##############################################
######################## MLE for F distribution #####################

n=500
y<-rf(n,10,5)
n <- length(y)
MLE.f<-function(df0,y){  
  floglike <- function(df,y) {
    f<- outer(y,df,FUN= function(y,df) df(y,df,5,log = TRUE))
    f<-(-1)*apply(f,2,sum)
    return(f)
  }
  result<-optim(df0,floglike,y=y,method = "Brent", lower = 0, upper = 50, hessian = TRUE)
  return(result)
}

Re<-MLE.f(2,y)
Re

df.hat<-Re$par
df.hat

# bootstraping for find distribution of MLE

B=2000   
n<-length(y)
dfhat.Bootstrap<-numeric(B)
for(i in 1:B){
  ystar<-rf(n,df.hat,5)
  dfhat.Bootstrap[i]<-MLE.f(0.5,ystar)$par
}

dfhat.Bootstrap


# histogram
hist(dfhat.Bootstrap, main = "sampling distribution of Bootrapped MLE for fisher", xlab = "estimated value", breaks = 20, prob = TRUE)

abline(v=mean(dfhat.Bootstrap),col=290,lty=2,lwd=2)
abline(v=10,col="green",lty=2,lwd=2)

lines(density(dfhat.Bootstrap),col="blue",lty=3,lwd=2)

mu<-mean(dfhat.Bootstrap)
sigma<-sd(dfhat.Bootstrap)

curve(dnorm(x,mu,sigma),col=colors()[91],add=TRUE)


### Another Version of Bootstrap 
B=2000   
n<-length(y)
dfhat.Bootstrap1<-numeric(B)
for(i in 1:B){
  ystar1<-sample(y,n,replace=TRUE)
  dfhat.Bootstrap1[i]<-MLE.f(2,ystar1)$par
}


# histogram
hist(dfhat.Bootstrap1, main = "sampling distribution of Bootrapped MLE for fisher", xlab = "estimated value", breaks = 20, prob = TRUE,col="gray95")
abline(v=mean(dfhat.Bootstrap1),col=2,lty=2,lwd=2)
lines(density(dfhat.Bootstrap1),col="blue",lty=3)
mu<-mean(dfhat.Bootstrap1)
sigma<-sd(dfhat.Bootstrap1)
curve(dnorm(x,mu,sigma),col=colors()[91],add=TRUE)

### comparing methods 
plot(density(dfhat.Bootstrap),ylim=c(0,0.4),col="blue",lty=2)
lines(density(dfhat.Bootstrap1),col="red",lty=3)
abline(v=mean(dfhat.Bootstrap),col= "blue",lty=2)
abline(v=mean(dfhat.Bootstrap1),col= "red",lty=3)
abline(v=10,col= "green",lty=3)



par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(dfhat.Bootstrap,col="gray70")
abline(h=mean(dfhat.Bootstrap),col="red",lty=2)

ts.plot(dfhat.Bootstrap1,col="gray70")
abline(h=mean(dfhat.Bootstrap1),col="red",lty=2)


# confidence interval by bootstraping methods

# Hessian methods
Re<-MLE.f(2,y)
Re
df.hat<-Re$par
SE.df<-sqrt((Re$hessian)^(-1))
alpha=0.05
L<-df.hat-qnorm(1-alpha/2)*SE.df
L
U<-df.hat+qnorm(1-alpha/2)*SE.df
U
c(L,U)

# Bootstrap methods

conf.boot<-quantile(dfhat.Bootstrap, probs=c(0.025,0.975), names=FALSE) 
conf.boot
conf.boot1<-quantile(dfhat.Bootstrap1, probs=c(0.025,0.975), names=FALSE) 
conf.boot1

## plot confidence intervals

par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(dfhat.Bootstrap,col="gray70")
abline(h=mean(dfhat.Bootstrap),col="red",lty=2)
abline(h=c(L,U),lty=3,col="blue")



ts.plot(dfhat.Bootstrap1,col="gray70")
abline(h=mean(dfhat.Bootstrap1),col="red",lty=2)
abline(h=conf.boot,lty=3,col="blue")
abline(h=conf.boot1,lty=3,col="green")


polygon(c(1:B, rev(1:B)), c(rep(conf.boot[1],B), rev(rep(conf.boot[2],B))), col = rgb(.9,.2,.8,alpha=.2), border = NA)

########################################################
######################################
##The Log Chi-Squared Distribution##
#"df" as parameter#
set.seed(143)
n=100
y<-rchisq(n,2) # my orginal "y"
Nloglike.ch<-function(df,y){
  f<-outer(y,df,FUN=function(y,df) dchisq(y,df,log=TRUE))
  f<-(-1)*apply(f,2,sum) # "-1" made my log likelihood negative
  return(f)
}
df<-1:3
Nloglike.ch(df,y)
y
hist(y)
result<-optim(0.5, Nloglike.ch,y=y,method="Brent",lower = 0.05,upper = 10,hessian = TRUE)
result
df.hat<-result$par
df.hat #show parameter
se.mle<-sqrt(result$hessian^(-1))
se.mle
# frome library(knitr)
print(kable(cbind(direct=c(mean=mean(y), sd=sqrt(mean(y))), optim =c(result$par,se.mle)),digits = 3))


# MLE function for Chi-Squared Distribution 
set.seed(143)
n=100
y<-rchisq(n,2) # my orginal "y"
MLE.chisq<-function(df0,y){  
  Nloglike.ch<-function(df,y){
    f<-outer(y,df,FUN=function(y,df) dchisq(y,df,log=TRUE))
    f<-(-1)*apply(f,2,sum) # "-1" made my log likelihood negative
    return(f)
  }
  result<-optim(df0,Nloglike.ch,y=y,method = "Brent", lower = 0.1, upper = 10, hessian = TRUE)
  return(result)
}

Re<-MLE.chisq(0.5,y)
Re
df.hat<-Re$par
df.hat

B=2000   ## iteration Bootstrap
n<-length(y)
Lhat.Bootstrap<-numeric(B)
for(i in 1:B){
  ystar<-rchisq(n,df.hat)
  Lhat.Bootstrap[i]<-MLE.chisq(0.5,ystar)$par
}
# histogram
hist(Lhat.Bootstrap, main = "sampling distribution of Bootrapped MLE for chisq", xlab = "estimated value", breaks = 10, prob = TRUE,col="gray93")

abline(v=mean(Lhat.Bootstrap),col="goldenrod4",lty=2,lwd=2)

lines(density(Lhat.Bootstrap),col="maroon",lty=2,lwd=2)

abline(v=2,col="green",lty=2,lwd=2)
mu<-mean(Lhat.Bootstrap)
sigma<-sd(Lhat.Bootstrap)

curve(dnorm(x,mu,sigma),col="gray0",add=TRUE)


### Another Version of Bootstrap 
B=2000   ## iteration Bootstrap
n<-length(y)
Lhat.Bootstrap1<-numeric(B)
for(i in 1:B){
  ystar1<-sample(y,n,replace=TRUE)
  Lhat.Bootstrap1[i]<-MLE.chisq(0.5,ystar1)$par
}


# histogram
hist(Lhat.Bootstrap1, main = "sampling distribution of Bootrapped MLE for normal", xlab = "estimated value", breaks = 20, prob = TRUE,col="gray95")
abline(v=mean(Lhat.Bootstrap1),col=2,lty=2,lwd=2)
lines(density(Lhat.Bootstrap1),col="blue",lty=3)
mu<-mean(Lhat.Bootstrap1)
sigma<-sd(Lhat.Bootstrap1)
curve(dnorm(x,mu,sigma),col=colors()[91],add=TRUE)

### comparing methods 
plot(density(Lhat.Bootstrap),col="blue",lty=2)
lines(density(Lhat.Bootstrap1),col="red",lty=3)
abline(v=mean(Lhat.Bootstrap),col= "blue",lty=2)
abline(v=mean(Lhat.Bootstrap1),col= "red",lty=3)

par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(Lhat.Bootstrap,col="gray70")
abline(h=mean(Lhat.Bootstrap),col="red",lty=2)

ts.plot(Lhat.Bootstrap1,col="gray70")
abline(h=mean(Lhat.Bootstrap1),col="red",lty=2)

# confidence interval by bootstraping methods

# Hessian methods
Re<-MLE.chisq(0.5,y)
Re
df.hat<-Re$par
SE.df<-sqrt((Re$hessian)^(-1))
alpha=0.05
L<-df.hat-qnorm(1-alpha/2)*SE.df
L
U<-df.hat+qnorm(1-alpha/2)*SE.df
U
c(L,U)

# Bootstrap methods

conf.boot<-quantile(Lhat.Bootstrap, probs=c(0.025,0.975), names=FALSE) 
conf.boot
conf.boot1<-quantile(Lhat.Bootstrap1, probs=c(0.025,0.975), names=FALSE) 
conf.boot1

## plot confidence intervals

par(mfrow=c(2,1),mar=c(2,2,2,2))
ts.plot(Lhat.Bootstrap,col="gray70")
abline(h=mean(Lhat.Bootstrap),col="red",lty=2)
abline(h=c(L,U),lty=3,col="orange")

ts.plot(Lhat.Bootstrap1,col="gray70")
abline(h=mean(Lhat.Bootstrap1),col="red",lty=2)
abline(h=conf.boot,lty=3,col="orange")
abline(h=conf.boot1,lty=3,col="green")

polygon(c(1:B, rev(1:B)), c(rep(conf.boot[1],B), rev(rep(conf.boot[2],B))), col = rgb(.9,.2,.8,alpha=.2), border = NA)

#################Beta logliklihood with two parameter#################
n=100
set.seed(143)
y<-rbeta(n,2,5)
MLE.Beta<-function(alpha0,Beta0,y){
Loglik.Beta<-function(alpha,Beta,y){
  g<-0
  for(i in 1:length(alpha)){
    f<-outer(Beta,y,FUN=function(Beta,y) dbeta(y,alpha[i],Beta,log=TRUE))
    f<-(-1)*apply(f,1,sum)
    g<-c(g,f)
   
  }
  g<-g[-1]
  return(g)
}
result<-optim(c(alpha0,Beta0),function(x) Loglik.Beta(alpha=x[1],Beta=x[2],y),hessian = TRUE)
return(result)

}


Re<-MLE.Beta(6,7,y)
Re

alpha.hat<-Re$par[1]
alpha.hat

Beta.hat<-Re$par[2]
Beta.hat


S<-solve(Re$hessian)
S
SE.alpha<-sqrt(diag(S)[1])
SE.alpha
alpha=0.05
L<-alpha.hat-qnorm(1-alpha/2)*SE.alpha
L
U<-alpha.hat+qnorm(1-alpha/2)*SE.alpha
U
c(L,U)

SE.Beta<-sqrt(diag(S)[2])
SE.Beta
L1<-Beta.hat-qnorm(1-alpha/2)*SE.Beta
L1
U1<-Beta.hat+qnorm(1-alpha/2)*SE.Beta
U1
c(L1,U1)
###################################################
#### class 18 #### # simulation book
###################################################
##example 1-1 , page"1"##
sample(1:100,5,replace = F)
x<-sample(1:100,5,replace = F)
x
x<=90
sum(x<=90)

iter<-1000000
N<-0
for (i in 1:iter) {
  x<-sample(1:100,5,replace = F)
  N[i]<-sum(x<=90)
}
#N
p<-mean(N==5)
p
hist(N)
ts.plot(N)

Con<-function(iter){
  N<-numeric(iter)
  for (i in 1:iter) {
    x<-sample(1:100,5,replace = F)
    N[i]<-sum(x<=90)
  } 
  p<-mean(N==5)
  return(p)
}
Con(10)
Con(100)
Con(1000)
Con(100000)
n<-seq(100,10^6,85000)
n
tic1 <- proc.time()
Prob<-0
for (i in 1:length(n))  {
  Prob[i]<-Con(n[i])
}
Prob
toc1<-proc.time()
toc1[3]-tic1[3]
plot(n,Prob,type = "l")
points(n,Prob,pch=16,col=2)
abline(h=0.5838,lty=2,col=4)

## with lapply instead for
n<-seq(100,10^6,85000)
n
tic <- proc.time()
FUN1<-function(i) sum(sample(1:100,5,rep=F)<=90)
FUN2<-function(iter) mean(as.numeric(lapply(1:iter,FUN1))==5)
#FUN2(1000)
Prob<-as.numeric(lapply(n,FUN2))
Prob
toc<-proc.time()
toc[3]-tic[3]
plot(n,Prob,type = "l")
points(n,Prob,pch=16,col=2)
abline(h=0.5838,lty=2,col=4)

####### with one for and lapply
tic2 <- proc.time()
Prob1<-as.numeric(lapply(n,Con))
Prob1
toc2<-proc.time()
toc2[3]-tic2[3]
##################################
#example 2: toss dice (see six number? what is prob?)
# ,p(x=6)=?  #
real=1/6
sample(1:6,1)
x<-sample(1:6,1)
x<-0
for(i in 1:1000000){
  x[i]<-sample(1:6,1)
  
}
#x
p6<-mean(x==6)
p6
real
real-p6

sample(1:6,10,replace = T)
y<-sample(1:6,1000000,replace = T)
p6<-mean(x==6)
p6



n<-seq(10,10^6,1000)
n
Prob<-0
for (i in 1:length(n))  {
  Prob[i]<-mean(sample(1:6,n[i],replace = T)==6)
}
Prob
plot(n,Prob,type = "l")
points(n,Prob,pch=16,col=2)
abline(h=real,lty=2,col=4)

## problem 1:coin tossing, or heads or tails?
## problem 2: Coverage Probabilities of Binomial Confidence Intervals 

#############################################
##### class 19 #####
#############################################
##### example1-2 page:4 birthday maching ####
x<-sample(1:365,25,replace = TRUE)
x
y<-unique(x)
length(x)
length(y)
length(x)-length(y)
table(x)
BM<-function(iter){
  n<-numeric(iter)
  for(i in 1:iter){
  x<-sample(1:365,25,replace = TRUE)
  y<-unique(x)
  n[i]<-length(x)-length(y)
  }
  return(n)
}
w<-BM(10000)
#mean(w==0)
p<-1-mean(w==0)
p
real<-0.5687
real-p
#### convergence plot ####
iter<-seq(1000,10^6,20000)
#iter
for(i in 1:length(iter)){
  w<-BM(iter[i])
  p[i]<-1-mean(w==0)
}
p
plot(iter,p,type="l",ylim = c(0.561,0.579))
points(iter,p,pch=16,cex=0.3,col=2)
abline(h=real,v=400000,lty=2,col=4)
cbind(iter,p)
phat.converge<-mean(p[-(1:20)])
abline(h=phat.converge,lty=3,col=2)
phat.converge
real
round(phat.converge,4)

#############################################
#######  class 20 #################
##############################################
#generating random numbers
#exampel 2.1 ,page 25
d<-53
b<-0
a<-20
n<-20
r<-numeric(n)
#r
seed=12
r[1]<-seed
for(i in 1:(n-1)) r[i+1]=(a*r[i]+b)%%d 
r
u<-(r+0.5)/d
u
hist(u)

u1<-u[1:(n-1)]
u2<-u[2:n]
plot(u1,u2,pch=20) 

#example 2.2,page28 (IBM)
d<-86436
b<-18257
a<-1093
n<-1000
r<-numeric(n)
seed=7
r[1]<-seed
for(i in 1:n){
r[i+1]=(a*r[i]+b)%%d
}
u<-(r+0.5)/d
hist(u)
u1<-u[1:(n-1)]
u2<-u[2:n]
plot(u1,u2,pch=20) 

plot(u1,u2,xlim=c(0,0.1),ylim=c(0,0.1),pch=20)

#example2.3
d<-2^31
b<-0
a<-65539
n<-5000 
r<-numeric(n)
seed=7
r[1]<-seed
for(i in 1:n){ r[i+1]=(a*r[i]+b)%%d }
u<-(r+0.5)/d
hist(u)
u1<-u[1:(n-1)]
u2<-u[2:n]
plot(u1,u2,pch=20) 
plot(u1,u2,xlim=c(0,0.02),ylim=c(0,0.02),pch=20)
#######

###########################
# X=Y in distribution

#use Q-Q plot
x<-rnorm(1000)
y<-rt(1000,2)
#H0: X=Y in distribution
qqplot(x,y)
abline(0,1,col=2,lty=2) 
curve(dnorm(x),-10,10,ylim=c(0,0.6)) 
curve(dt(x,2),-10,10,col=2,add=TRUE)

# ks.test
ks.test(x,y)
###################################

qqplot(u,runif(2000),pch=20,cex=.1)

ks.test(u,"punif",0,1)
###################################### 
#x~exp(lambda)
#generating random numbers of x~exp(lambda) # E(X)=lambda
lambda=2
n=1000
u<-runif(n)
u 
x<-(-lambda*log(1-u))
 x

mean(x) 
hist(x,freq=FALSE)
abline(v=mean(x),col=2,lty=3)
abline(v=lambda,col=5,lty=3) 
#text(10,0.3,"mean of x",col=25,cex=2) 
legend("topright",c("hist of x","mean of x","lambda"),col=c(1,2,5),lty=c(1,3,3))
lines(density(x),col=2,lty=2) 
curve(1/lambda*exp(-x/lambda),0,20,col=6,cex=2,add=TRUE) #############################
ppoints(10)
n=100000
u<-runif(n)
x<-(-lambda*log(1-u))
#x
qqplot(qexp(ppoints(100),2),x) 
qqline(x,distribution=function(p) qexp(p,2))
## t-student?
qqplot(qt(ppoints(100),4),x)
qqline(x,distribution=function(p) qt(p,4)) ############################################

############################################ 
##x~Gamma(alpha,1) #y^(alpha-1)*exp(-y)/gamma(alpha)
#alpha=10, 5, 2 

#############################################
#######  class 21 #################
##############################################
## Exercise: generate samples from the following distributions
# by CDF method

#1: laplace distribution f(x)=1/2*exp(-|x|)
#2:  chi-square distribution df=2 and 4


#Example: Use CDF method to simulate a random variable which pdf is f(x)=3*x^2, 0<x<1




# Example:
#The cdf is given by F(x)=x^3
#so, the inverse is x=F^−1(u)=u^(1/3)

# Checking if it went well
n=10000
u<-runif(n)
vals<-u^(1/3)
hist(vals, breaks=50, freq=FALSE, main=expression("Sample vs true Density [ f(x)=" ~3*x^2~"]"))
curve(3*x^2, 0, 1, col="red", lwd=2, add=T)


######################################## 
#box_muller transformation 
#########################################
#to generat random sample from normal distribution #########################################
n=10000
u1<-runif(n)
u2<-runif(n)
z1<-sqrt(-2*log(u1))*cos(2*pi*u2)
#z1
z2<-sqrt(-2*log(u1))*sin(2*pi*u2)
#z2
par(mfrow=c(2,2))
hist(z1,freq=FALSE)
lines(density(z1),col=2,lty=2) 
curve(dnorm(x),-4,4,col=3,lty=3,add=TRUE) 
qqnorm(z1)
qqline(z1)
hist(z2,freq=FALSE)
lines(density(z2),col=2,lty=2) 
curve(dnorm(x),-4,4,col=3,lty=3,add=TRUE) 
qqnorm(z2)
qqline(z2)
shapiro.test(z1)
shapiro.test(z2)

########################################################
#Exercise
##Z=N(0,1) ===>  x=Z^2~?  ,  y=sqrt(x)~?

# Transformation of distributions
#########################################################
#chapter 3
########################
# Riemann approximation:
########################
#integrate e^-x in [0,1]
a=0
b=1
n=1000
w<-(b-a)/n
ksi<-seq(a+w/2,b-w/2,length=n)
h<-exp(-ksi)
s<-sum(h*w)
s
real<-1-exp(-1)
real
error<-real-s
error
#####integrate e^(-x^2) in [0,1]
a=0
b=1
n=1000
w<-(b-a)/n
ksi<-seq(a+w/2,b-w/2,length=n)
h<-exp(-ksi^2)
s<-sum(h*w)
s
real<-integrate(function(x) exp(-x^2),a,b)
real
#names(real)
error<-real$value-s
error
#?integrate
################################
# Monte Carlo approximation
################################ 
 #Monte Carlo integration is an estimation of the true integration based on random sampling and in the Strong Law of Large Numbers.

#The estimator of θ=∫g(x)dx on [a,b] is computed as follows:

#1:Generate X1,X2,…,Xn iid from Unif(a,b)
#2:Compute g_hat=mean(g(X1,..,Xn))
#The estimation θ̂=(b−a)g_hat
#θ̂ is itself a random variable, which by the Strong Law of Large Numbers: θ̂ →θ as n→∞
##################################

# Monte Carlo Method:
# p(0<z<1)=?  z~N(0,1)
n=1000
a=0
b=1
u<-runif(n)
ghat<-mean(dnorm(u))
(b-a)*ghat


#random sampleing 
## z ~ N(0, 1)
# p(0<z<1) = ?   is ∫f(x)dx on [0,1]
n=100000
z<-rnorm(n)
p<-mean((0<=z) & (z<=1))
p

## Example: f(x)=e^(-x^2) in [0,1]
#E(f(U))=?
f<-function(x) exp(-x^2)
a=0
b=1
n=1000
u<-runif(n,a,b)

(b-a)*mean(f(u))

## Example: f(x)=e^(-x^2) in [-1,0]
#E(f(U))=?
f<-function(x) exp(-x^2)
a=-1
b=0
n=10000
u<-runif(n,a,b)

(b-a)*mean(f(u))
curve(f,-1,1)
abline(v=c(-1,0),lty=2,col=2)



## Exercise
#1: x ~ exp(lambda = 5)
#2:x ~ chisq(4)
#3: x ~ t-student(5)
# p(1<x<5)
# p(10<x<100) , ...

# Riemann methods and MC?

###############################################
#Exercise
#intrgerate cos(2*pi*x)*exp(-x*sin(2*pi*x))  -1<x<1
# Riemann and Monte Carlo approximations?

# Riemann approximation
a=-1
b=1
n=100
w=(b-a)/n
ksi<-seq(a+w/2,b-w/2,length=n)
h<-cos(2*pi*ksi)*exp(-ksi*sin(2*pi*ksi))
s<-sum(w*h)
s
real<-integrate(function(x) cos(2*pi*x)*exp(-x*sin(2*pi*x)),a,b) 
real
error<-real$value-s
error

# Monte Carlo approximation ?

##########################################################
#generating sample random from f(x)=1/2*exp(-|x|) by CDF
n=10000
u<-runif(n)
#?ifelse
x<-ifelse(u<=0.5,log(2*u),-log(2*(1-u)))
 hist(x,freq=FALSE,ylim=c(0,0.5))
lines(density(x),lty=2,col=2)
 curve(1/2*exp(-abs(x)),-6,6,col=4,add=TRUE) ###############################################


###############################
## Acceptance-Rejection Method
###############################
# Reference 
#Robert, C.P., Casella, G., Introducing Monte Carlo Methods with R, New York: Springer (2010).

#Problem: Generate X∼f
# from an arbitray pdf f (especially when it’s hard to sample from f).

#We must find Y∼g under the only restriction that 
#∀f(x)>0:f(x)<c.g(x),c>1.

#Instead of sampling from f(x)
# which might be difficult, we use c.g(x)
#to sample instead.

#For each value required the method follows:

#Generate a random y from g
#Generate a random u from U(0,1)
#If u<f(y)/(c.g(y)) then return y
#else reject y
# and goto 1.

#The algorihtm in R:

# generate n samples from f using rejection sampling with g (rg samples from g)
accept.reject <- function(f, c, g, rg, n) { 
  n.accepts     <- 0
  result.sample <- rep(NA, n)
  
  while (n.accepts < n) {
    y <- rg(1)               # step 1
    u <- runif(1,0,1)        # step 2
    if (u < f(y)/(c*g(y))) { # step 3 (accept)
      n.accepts <- n.accepts+1
      result.sample[n.accepts] = y
    }
  }
  
  result.sample
}
###########################################
#Example: generate samples from distribution Beta(2,2), where we use the uniform has g since f(x)<2*g(x):

f  <- function(x) 6*x*(1-x)     # pdf of Beta(2,2), maximum density is 1.5
g  <- function(x) x/x           # g(x) = 1 but in vectorized version
rg <- function(n) runif(n,0,1)  # uniform, in this case
c  <- 2                         # c=2 since f(x) <= 2 g(x)

vals <- accept.reject(f, c, g, rg, 10000) 

# Checking if it went well
hist(vals, breaks=30, freq=FALSE, main="Sample vs true Density")
xs <- seq(0, 1, len=100)
lines(xs, dbeta(xs,2,2), col="red", lwd=2)



##########################
#### package AR
##########################
# install.packages("AR")
library(AR)

# Example 1:f(x)=beta(2,2) ? 0<x<1
data = AR.Sim( n = 100, f_X = function(y){dbeta(y,2,2)}, Y.dist = "unif", Y.dist.par = c(0,1), Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# QQ-plot
q <- qbeta(ppoints(100), 2, 2)
qqplot(q, data, cex=0.6, xlab="Quantiles of Beta(2,2)",
       ylab="Empirical Quantiles of simulated data")
abline(0, 1, col=2)


# Example 2:f(x)=beta(2,2) ? 0<x<1,  g(x)=U(0,0.5)
data1 = AR.Sim( n = 100, f_X = function(y){dbeta(y,2,2)}, S_X=c(0,1),Y.dist = "unif", Y.dist.par = c(0,0.5), Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# QQ-plot
q <- qbeta(ppoints(100), 2, 2)
qqplot(q, data1, cex=0.6, xlab="Quantiles of Beta(2,2)",
       ylab="Empirical Quantiles of simulated data")
abline(0, 1, col=2)

# Example 3:f(x)=beta(2,2) ? 0<x<1   g(x)=N(0.5,0.3)
data2= AR.Sim( n = 100, f_X = function(y){dbeta(y,2,2)}, Y.dist = "norm", Y.dist.par = c(.5,.3), Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# QQ-plot
q <- qbeta(ppoints(100), 2, 2)
qqplot(q, data2, cex=0.6, xlab="Quantiles of Beta(2,2)",
       ylab="Empirical Quantiles of simulated data")
abline(0, 1, col=2)


## 

#Example2: f(x)=1/2 exp(-|x|)   -inf<x<inf
f<-function(x) 1/2*exp(-abs(x))

#g?

curve(f,-6,6,col=4) 
curve(dnorm,-6,6,col=2,add=T)
curve(dt(x,3),-6,6,col=6,add=T)

#g=U(-6,6)
data1 = AR.Sim( n = 100, f_X = function(x) {1/2*exp(-abs(x))}, Y.dist = "unif", Y.dist.par = c(-6,6), xlim=c(-6,6),Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# Checking if it went well
hist(data1, breaks=30, freq=FALSE, main="Sample vs true Density")
curve(f,-6,6,col="Orange1",add=TRUE) 


#g=U(-10,10)
data2 = AR.Sim( n = 100, f_X = function(x) {1/2*exp(-abs(x))}, Y.dist = "unif", Y.dist.par = c(-10,10), xlim=c(-10,10),Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# Checking if it went well
hist(data2, breaks=30, freq=FALSE, main="Sample vs true Density")
curve(f,-6,6,col="Orange1",add=TRUE) 


#g=N(0,10)
data3 = AR.Sim( n = 100, f_X = function(x) {1/2*exp(-abs(x))}, Y.dist = "norm", Y.dist.par = c(0,10), xlim=c(-6,6),Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# Checking if it went well
hist(data3, breaks=30, freq=FALSE, main="Sample vs true Density")
curve(f,-10,10,col="Orange1",add=TRUE) 



#g=t(3)
data4 = AR.Sim( n = 1000, f_X = function(x) {1/2*exp(-abs(x))}, Y.dist = "t", Y.dist.par = 3, xlim=c(-10,10),Rej.Num = TRUE, Rej.Rate = TRUE, Acc.Rate = FALSE )
# Checking if it went well
hist(data4, breaks=30, freq=FALSE, main="Sample vs true Density")
curve(f,-10,10,col="Orange1",add=TRUE) 

######################################################################
#######  class 22 #################
#####################################################################




# Example:given two standard normal iid random vars X1,X2 estimate E[|X1−X2|]
# Monte Carlo
n <- 1e5  # is 10000
X1 <- rnorm(n,0,1)
X2 <- rnorm(n,0,1)

thetas <-abs(X1-X2)
mean(thetas)

mean(abs(X1-X2))

mean(X1)
mean(X2)
var(X1)
# standard error = sqrt( var(|X1-X2|)/n )
sqrt(var(thetas)/n)