1

I have a problem regarding the following model,

enter image description here

where I want to make inference on μ and tau, u is a known vector and x is the data vector. The log-likelihood is

enter image description here

I have a problem writing a log-likelihood in R.

x <- c(3.3569,1.9247,3.6156,1.8446,2.2196,6.8194,2.0820,4.1293,0.3609,2.6197)
mu <- seq(0,10,length=1000)

normal.lik1<-function(theta,x){ 
  u <- c(1,3,0.5,0.2,2,1.7,0.4,1.2,1.1,0.7)  
  mu<-theta[1] 
  tau<-theta[2] 
  n<-length(x) 

logl <-  sapply(c(mu,tau),function(mu,tau){logl<- -0.5*n*log(2*pi) -0.5*n*log(tau^2+u^2)- (1/(2*tau^2+u^2))*sum((x-mu)^2) } )

  return(logl) 
  } 

#test if it works for mu=1, tau=2
head(normal.lik1(c(1,2),x))
#Does not work..

I want to be able to plug in the vector for mu and plot it over mu for a fixed value of tau, say 2. I also want to find out the MLE's of tau and mu using the optim function. I tried:

theta.hat<-optim(c(1,1),loglike2,control=list(fnscale=-1),x=x,,method="BFGS")$par

But it does not work.. Any suggestions to how I can write the likelihood?

Hans Christensen
  • 39
  • 8

1 Answers1

2

First, as has been mentioned in the comments to your question, there is no need to use sapply(). You can simply use sum() – just as in the formula of the logLikelihood.

I changed this part in normal.lik1() and multiplied the expression that is assigned to logl by minus 1 such that the function computes the minus logLikelihood. You want to search for the minimum over theta since the function returns positive values.

x < c(3.3569,1.9247,3.6156,1.8446,2.2196,6.8194,2.0820,4.1293,0.3609,2.6197)
u <- c(1,3,0.5,0.2,2,1.7,0.4,1.2,1.1,0.7) 

normal.lik1 <- function(theta,x,u){ 
  mu <- theta[1] 
  tau <- theta[2] 
  n <- length(x) 
  logl <- - n/2 * log(2*pi) - 1/2 * sum(log(tau^2+u^2)) - 1/2 * sum((x-mu)^2/(tau^2+u^2))
  return(-logl) 
}

This can be done using nlm(), for example

nlm(normal.lik1, c(0,1), hessian=TRUE, x=x,u=u)$estimate

where c(0,1) are the starting values for the algorithm.

To plot the logLikelihood for a range of values of mu and some fixed tau you can adjust the function such that mu and tau are separate numeric arguments. 

normal.lik2 <- function(mu,tau,x,u){ 
  n <- length(x) 
  logl <- - n/2 * log(2*pi) - 1/2 * sum(log(tau^2+u^2)) - 1/2 * sum((x-mu)^2/(tau^2+u^2))
  return(logl) 
}

Then define some range for mu, compute the loglikelihood and use plot().

range.mu <- seq(-10,20,0.1)

loglik <- sapply(range.mu, function(m) normal.lik2(mu=m,tau=2,x=x,u=u))

plot(range.mu, loglik, type = "l")

enter image description here

I'm sure there are more elegant ways to do this but this does the trick.

Martin C. Arnold
  • 9,483
  • 1
  • 14
  • 22
  • makes sense! If I were to change the -logL to logL in normal.lik1, would then: "optim(c(1,1),normal.lik1,control=list(fnscale=-1),x=x, u=u,method="BFGS")$par" do the exact same thing as nlm? – Hans Christensen Jan 10 '18 at 12:22
  • Admittedly I'm not too familiar with `optim()` and `nlm()` let alone the algorithms used but for the data at hand both yield estimates that are virtually the same (mu = 2.1863, tau = 0). – Martin C. Arnold Jan 10 '18 at 14:34
  • I know, whats bothers me is that I dont think tau should be 0, in fact I assume it should be between 0.5 and 3. The initial model is X_i=mu+lambda_i+epsilon_i. where lambda is N(0,tau^2) and epsilon_i is N(0,u^2_i). Therefore I assumed X is N(mu,lambda^2+u^2). That means tau cant be negative, and I would assume the MLE for it should take on a value bewteen 0 and say 5... perhaps I got the likelihood wrong? – Hans Christensen Jan 10 '18 at 15:07
  • I get the same mu and tau from using optim BTW – Hans Christensen Jan 10 '18 at 15:08
  • Your derivation of the logLik is correct. However, there was another error in normal.lik1() (second summand multiplied by `n`). I edited my answer. Now I get (mu=2.6525, tau=0.7692) using nlm(). – Martin C. Arnold Jan 10 '18 at 16:44