# assume X = c(X_1,...,X_n) an iid sample from f(x|theta)
# here we estimate the two key quantities used in the development of the CRLB
# estimated expected slope of the log f(X|theta) = E[d/dtheta log f(X|theta)] = 0
# estimated variance of the slope of the log f(x|theta) = 
#     Var(d/dtheta log f(X|theta)) = E[(d/dtheta log f(X|theta))^2] = 
#                                   -E[d^2/dtheta^2 log f(X|theta)] = I(theta)

# Example 1.
# X ~ Exp(theta)
# f(x|theta) = theta^-1 exp(-x/theta) 0 < x

# L(theta) = theta^-n exp(-sum(x_i)/theta)
# l(theta) = log L(theta) = -n log(theta) - sum(x_i) theta^-1
# d.l(theta) = -n theta^-1 + sum(x_i) theta^-2
# d2.l(theta) = n theta^-2 - 2 sum(x_i) theta^-3

# I(theta) = n theta^-2
# AV = theta^2/n

# simulate some date

theta.0 = 3; lambda.0 = 1/theta.0
n = 100

x = rexp(n,rate = lambda.0)
x

# plot the fitted model on the histogram

hist(x,freq=FALSE)

theta.mle = mean(x)
lambda.mle = 1/theta.mle

curve(dexp(x,lambda.mle),add=T)

# plot the functions of theta

theta = seq(0.01,6,0.1)

par(mfrow=c(2,2))

# Likelihood

Lik = function(theta){
  theta^-n * exp(-sum(x)/theta)
}

plot(theta,Lik(theta),type='l')

# log Likelihood

log.Lik = function(theta){
  -n*log(theta) - sum(x)/theta
}

plot(theta,log.Lik(theta),type='l')

# derivative of the log.Lik

d.log.Lik = function(theta){
  -n/theta + sum(x)/(theta^2)
}

plot(theta,d.log.Lik(theta),type='l')

# second derivative of the log.Lik

d2.log.Lik = function(theta){
  n/(theta^2) - 2*sum(x)/(theta^3)
}

plot(theta,d2.log.Lik(theta),type='l')

# estimates

mean(d.log.Lik(theta))

var(d.log.Lik(theta))