%auto-ignore #### function to compute log likelihood for MVN NLL = function(MU,SIG,XY){ SIGINV = solve(SIG) ll = 0 for(ii in 1:dim(XY)[1]){ ll = ll + t(matrix(XY[ii,],ncol=1) - MU)%*%SIGINV%*% (matrix(XY[ii,],ncol=1) - MU) } return(as.numeric(-ll/2)) } #### function to compute test statistic for MVN #### note we could have used -2*( NLL(MU,SIGHAT,XY) - NLL(MU0,SIGHAT,XY)) #### instead of N*t(MU - MU0)%*%SIGINV%*%(MU - MU0) T = function(MU0,MU,SIGINV,N){ ll = N*t(MU - MU0)%*%SIGINV%*%(MU - MU0) return(as.numeric(ll)) } ##### function to compute log likelihood ratio minus critical value Tsoln = function(d,phi,MU,SIGINV,N,chi2){ MU0 = MU + matrix(c(d*cos(phi),d*sin(phi)),ncol=1) return(T(MU0,MU,SIGINV,N) - chi2) } #### necessary library to compute bivariate normal data library(MASS) #### seed set.seed(2004) #### example 1 from Section 3.3 ##### data generation XY = mvrnorm(10,c(0,0),matrix(10*c(1,1/2,1/2,1),nrow=2)) #### plots data plot(XY,xlab='',ylab='',pch=20) ### maximum likelihood estimators MUHAT = matrix(c(mean(XY[,1]),mean(XY[,2])),ncol=1) SIGHAT = var(XY) SIGHATINV = solve(SIGHAT) ### function value at MLEs fmle0 = NLL(MUHAT,SIGHAT,XY) #### this holds radial coordinate r for each phi r_phi = numeric(180) #### this holds cartesian coordinate for each phi ci = matrix(0,nrow=180,ncol=2) ### set of angles to compute bound ### will ignore 181st phi = seq(0,2*pi,length=181) ### saves time of each computation t_phi = numeric(180) for(i in 1:180){ ### starts computation time a = as.numeric(Sys.time()) ### finds radial angle g = uniroot(Tsoln,c(0,3),phi = phi[i],MU=MUHAT,SIGINV=SIGHATINV, N=N,chi2 = qchisq(.9,2)) ### saves computation time t_phi[i] = (as.numeric(Sys.time()) - a) ### save radial coordinate & cartesian coordinates r_phi[i] = g$root ci[i,] = MUHAT + r_phi[i]*c(cos(phi[i]),sin(phi[i])) } #### adds boundary to plot lines(rbind(ci,ci[1,])) ### grid method; computes ll over [-3,3]X[-3,3] ### x & y values xg = seq(-3,3,by=.05) yg = seq(-3,3,by=.05) ### saves log likelihood values ll = matrix(0,nrow=length(xg),ncol=length(yg)) ### saves computation time at each i,j entry t_grid = matrix(0,nrow=length(xg),ncol=length(yg)) for(i in 1:dim(ll)[1]){ for(j in 1:dim(ll)[2]){ ### extracts x,y point MU0 = matrix(c(xg[i],yg[j]),ncol=1) ### starts computation time a = as.numeric(Sys.time()) ### computes log likelihood ratio ll[i,j] = T(MU0,MUHAT,SIGHATINV,N) ### saves computation time t_grid[i,j] = (as.numeric(Sys.time()) - a) } } #### plots data plot(XY,xlab='',ylab='',pch=20) #### adds contour contour(xg,yg,ll,levels = qchisq(.9,2),drawlabels=FALSE,add=TRUE) ### computation times in seconds ### radial profile log likelihood sum(t_phi) ### grid method sum(t_grid) #### example 2 from Section 3.3 #### identical to example 1 except XY is multiplied by 10 XY = 10*XY plot(XY,xlab='',ylab='',pch=20) MUHAT = matrix(c(mean(XY[,1]),mean(XY[,2])),ncol=1) SIGHAT = var(XY) SIGHATINV = solve(SIGHAT) fmle0 = NLL(MUHAT,SIGHAT,XY) r_phi = numeric(180) ci = matrix(0,nrow=180,ncol=2) phi = seq(0,2*pi,length=181) t_phi = numeric(180) for(i in 1:180){ a = as.numeric(Sys.time()) g = uniroot(Tsoln,c(0,30),phi = phi[i],MU=MUHAT,SIGINV=SIGHATINV, N=N,chi2 = qchisq(.9,2)) t_phi[i] = (as.numeric(Sys.time()) - a) r_phi[i] = g$root ci[i,] = MUHAT + r_phi[i]*c(cos(phi[i]),sin(phi[i])) } lines(rbind(ci,ci[1,])) xg = seq(-30,30,by=.05) yg = seq(-30,30,by=.05) ll = matrix(0,nrow=length(xg),ncol=length(yg)) t_grid = matrix(0,nrow=length(xg),ncol=length(yg)) for(i in 1:dim(ll)[1]){ for(j in 1:dim(ll)[2]){ MU = matrix(c(xg[i],yg[j]),ncol=1) a = as.numeric(Sys.time()) ll[i,j] = -2*( NLL(MU,SIGHAT,XY) - fmle0) t_grid[i,j] = (as.numeric(Sys.time()) - a) } } plot(XY,xlab='',ylab='',pch=20) contour(xg,yg,ll,levels = qchisq(.9,2),drawlabels=FALSE,add=TRUE) sum(t_phi) sum(t_grid) ###### repeats example 2 100 times ###### set.seed(2525) nsim = 100 sigma = 1000 N = 10 phi = seq(0,2*pi,length=181) xg = seq(-30,30,by=.05) yg = seq(-30,30,by=.05) ### holds compuation times times = matrix(NA,nrow=nsim,ncol=2) ### holds boundaries using radial ciphi = array(NA,c(180,2,nsim)) ### holds all log likelihood ratio computation cigrd = array(NA,c(length(xg),length(yg),nsim)) for(z in 1:nsim){ XY = mvrnorm(N,c(0,0),matrix(sigma*c(1,1/2,1/2,1),nrow=2)) MUHAT = matrix(c(mean(XY[,1]),mean(XY[,2])),ncol=1) SIGHAT = var(XY) SIGHATINV = solve(SIGHAT) fmle0 = NLL(MUHAT,SIGHAT,XY) r_phi = numeric(180) ci = matrix(0,nrow=180,ncol=2) t_phi = numeric(180) for(i in 1:180){ a = as.numeric(Sys.time()) g = uniroot(Tsoln,c(0,75),phi = phi[i],MU=MUHAT,SIGINV=SIGHATINV, N=N,chi2 = qchisq(.9,2)) t_phi[i] = (as.numeric(Sys.time()) - a) r_phi[i] = g$root ci[i,] = MUHAT + r_phi[i]*c(cos(phi[i]),sin(phi[i])) } ciphi[,,z] = ci ll = matrix(0,nrow=length(xg),ncol=length(yg)) t_grid = matrix(0,nrow=length(xg),ncol=length(yg)) for(i in 1:dim(ll)[1]){ for(j in 1:dim(ll)[2]){ MU = matrix(c(xg[i],yg[j]),ncol=1) a = as.numeric(Sys.time()) ll[i,j] = -2*( NLL(MU,SIGHAT,XY) - fmle0) t_grid[i,j] = (as.numeric(Sys.time()) - a) } } cigrd[,,z] = ll times[z,] = c(sum(t_phi),sum(t_grid)) } #### radial profile log likelihood computation time sum(times[,1]) #### grid method computation time sum(times[,2])