# Edps/Psych/Stat 587
# HLM/Multilevel modeling
# Spring 2019
#
# Carolyn J. Anderson
#
# Lecture notes on Estimation
#
#

########################################################
# Small example of likelhood ratio:  n=1               #
########################################################

y <- 2
sigma2 <- 9
mu <- seq(from=-5, to=10, by=.1)

lr <- (1/sqrt(sigma2*pi*2)) * exp(-.5*(y-mu)**2/sigma2)

plot(mu,lr, type='l', col='blue',
     xlab = 'Possible Values for mu',
     ylab = 'Value of Likelihood Ratio Function',
     main = 'Likelihood Ratio Function:  y=2, sigma2=9')
          
       
 
########################################################
# A little bigger example of likelhood ratio:  n=5     #
########################################################
      
y <- c(-1, 2, 3, 6, 10)

lr5 <- c(1)
for (i in 1:5){
    lr5 = lr5* (1/sqrt(sigma2*pi*2))*exp(-.5*(y[i]-mu)**2/sigma2)
}
plot(mu,lr5, type='l', col='blue',
     xlab = 'Possible Values for mu',
     ylab = 'Value of Likelihood Ratio Function',
     main = 'Likelihood Ratio Function:  n=5, sigma2=9')
          
 
#######################################################
# Mu and Sigma2                                       #
#######################################################

sigma2 <- seq(from=.01, to=90, by=.1)

lr.bi <- matrix(1, nrow=151, ncol=900)


for (j in 1:151){
   for (k in 1:900){
     for (i in 1:5){
          lr.bi[j,k] = lr.bi[j,k]*(1/sqrt(sigma2[k]*pi*2))*exp(-.5*(y[i]-mu[j])**2/sigma2[k])
     }
  }
}


par(mfrow=c(1,1))
### 3 D plot ###

persp(mu,sigma2,lr.bi, col='lightpink', theta=45)


### countour ###
filled.contour(x = mu,
               y = sigma2,
               z = lr.bi,
               xlim = c(-2,9),
               ylim = c(0,60),
               xlab = 'Possible Values of Mu',
               ylab = 'Possible Value of Sigma2',
               main = 'Contour plot of likelihood ratio surface')


par(mfrow=c(1,2))
### 1 D plot of mu ####

one900 <- matrix(1, nrow=900, ncol=1)
mu.lr <- lr.bi %*% one900

plot(mu,mu.lr, type='l')


#### 1 D plot of sigma2 #######

one151 <- matrix(1, nrow=151,ncol=1)
s2.lr <- t(one151) %*% lr.bi

plot(sigma2,s2.lr, type='l', xlim= c(0, 100))