HW7comp

Ex 18_1:

There are 12 complicated crossroads in this city. The numbers of car accidents at the 12 crossroads within 1 month are as follows: 5, 5, 6, 0, 1, 2, 4, 4, 3, 5, 7, 10.

1. Please plot the log-likelihood function for mean number of car accidents within 1 month, lambda.

X <- c(5, 5, 6, 0, 1, 2, 4, 4, 3, 5, 7, 10)
lambda <- seq(0,10,0.001)
lambdahat <- mean(X)

loglik <- function(X,p){
        loglikelihood <- 0
        for(i in 1:length(X)){
                loglikelihood <-loglikelihood+ 
                        log((exp(-lambda)*(lambda^X[i]))/factorial(X[i]))
        }
        return(loglikelihood)
}

fp2 <- loglik(X,lambda)
plot(lambda,fp2,col=1, type = 'l', lwd = 2) 

abline(v = lambdahat, col = 2, lwd = 2)

2. Please solve the maximum likelihood estimate (MLE) of lambda by the Newton-Raphson method.And mark the MLE on the above plot.

decide the initial value by 0.8

lambda <- seq(0,10,0.001)
f <- (-length(X))+sum(X)/lambda
plot(lambda,f,type = 'l', lwd = 2)
abline(h = 0, col =2, lwd = 2)

define the function as below. the root 4.333 is the same from the mean of this dataset.

ftn2 <- function(lambda) {
        f <- (-length(X))+sum(X)/lambda
        #f <- (-10)+36/lambda # f = -n + summation of x*(1/lambda)
        df <- -(lambda^(-2))*sum(X) # df = -lambda^(-2) summation of x
        return(c(f, df))
}

newtonraphson(ftn2, 0.8, 1e-06)

## At iteration 1 value of x is: 1.452308 
## At iteration 2 value of x is: 2.417877 
## At iteration 3 value of x is: 3.486648 
## At iteration 4 value of x is: 4.1679 
## At iteration 5 value of x is: 4.327018 
## At iteration 6 value of x is: 4.333324 
## At iteration 7 value of x is: 4.333333 
## Algorithm converged

## [1] 4.333333

mean(X)

## [1] 4.333333

Ex S5_1 :

Data: BMIrepeated.csv on NTU COOL (Week 4). Please use “SEX”, “AGE”, “Treatment” as predictors, “BMI3” as response variable, to calculate the log likelihood of this model.

BMI <- read.csv("BMIrepeated.csv", stringsAsFactors = T)

# Constructing design matrix
X <- cbind(rep(1,length(BMI$BMI3)),BMI$SEX, BMI$Treatment, BMI$AGE)

# Calculating regression coefficients
beta <- solve(t(X)%*%X)%*%t(X)%*%matrix(BMI$BMI3,ncol = 1)

# Calculating residuals
re <- BMI$BMI3 - X%*%beta

# estimating sigma hat square
sigma2 <- sum(re^2)/length(BMI$BMI3)

-(length(BMI$BMI3)/2)*(log(2*pi)+log(sigma2)+1)

## [1] -337.3221

# linear model
model2 <- lm(BMI3 ~ SEX + AGE + Treatment, data = BMI)

# Calculating Log likelihood
logLik(model2)                                                          

## 'log Lik.' -337.3221 (df=5)