# Crops545 Lecture 17 Likelihood
# code provided by Zhiwu Zhang
# commented and modified by Yuanhong Song

# Univariate normal case ----

# X~N(mu, si)
# write a function: for given observation, erstimate the density (likelihood)
dnormal=function(x=0,mean=0,sd=1){
  p=1/(sqrt(2*pi)*sd)*exp(-(x-mean)^2/(2*sd^2))
  return(p)
}
# essentially, it does the same work as the R function dnorm()

# Say x=95 is observed, and X~ Norm(mu, si2)
# What could the value of mu and si2 possibly be?
x=c(95,95,95,95)
# here are some guess:
mean=c(100,100,85,85)
sd=c(1,2,5,10)

# use our own function to estimate the probability
dnormal(x,mean,sd)

# use the default R function to do the same thing
dnorm(x,mean,sd)   # R function

# Multiple variate normal case -----
# X~N-p(mu, SIGMA), given observation in p dimenssion, estimate likelihood
# the biggest difference here is the matrix calculation
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)
}

# create a multivariable data case
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)

va1=95
ve1=5
(V1=2*K*va1+ve1)
dmnormal(x,mean,V1)

va2=5
ve2=95
V2=2*K*va2+ve2
dmnormal(x,mean,V2)

va3=50
ve3=50
V3=2*K*va3+ve3
dmnormal(x,mean,V3)

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))


# The analysis above demonstrate the idea of maximum likelihood estimateGiven