Simulation to generate contour plot:
Two parameters: mean of g_1(t) = t, g_2(t)= I[t >= 0.5]. The true values are (1, exp(-0.5)).


library(emplik)
set.seed(123)
nn <- 1000
tvec <- rexp(nn)
cvec <- 0.1 + rexp(nn, rate=0.6)
yvec <- pmin(tvec, cvec)
dvec <- as.numeric(tvec <= cvec)
Cdata <- Wdataclean3(z=yvec, d=dvec)
mydd <- Cdata$dd
mydd[length(mydd)] <- 1

posi1 <- (mydd == 1)
#### KMcen <- WKM(x=Cdata$value, d=1-mydd)
KMcen1 <- WKM(x=Cdata$value, d=mydd)
NPmle <- sum((KMcen1$times)*(KMcen1$jump))
Alltime <- KMcen1$times[posi1]
Psdovec <- rep(0, nn)
Psdovec[posi1] <- Alltime*(nn*KMcen1$jump[posi1])
if (mean(Psdovec) != NPmle ) warning("somethingwrong")


mu1 <- 1:100/100 +0.6
mu2 <- 1:100/200 +0.4
myresult <- matrix(NA, ncol=100, nrow=100)
for(i in 1:100) for(j in 1:100) myresult[i,j]<-el.cen.EM2(x=yvec,d=dvec,fun=function(t){cbind(t,as.numeric(t>=0.5))},mu=c(mu1[i],mu2[j]))$"-2LLR"
Psdovec1 <- rep(0, nn)
Psdovec1[posi1] <- Alltime*(nn*KMcen1$jump[posi1])
Psdovec2 <- rep(0, nn)
Psdovec2[posi1] <- as.numeric(Alltime >= 0.5)*(nn*KMcen1$jump[posi1])
myresult2 <- matrix(NA, ncol=100, nrow=100)
for(i in 1:100)for(j in 1:100) myresult2[i,j]<-el.test(x=cbind(Psdovec1,Psdovec2),mu=c(mu1[i],mu2[j]))$"-2LLR"


contour(x=mu1, y=mu2, z=myresult2, level=c(1, 2, 3, 4, 5, 6, 7, 8), col="red")
par(new=TRUE)
contour(x=mu1, y=mu2, z=myresult, level=c(1, 2, 3, 4, 5, 6, 7, 8), col="blue")
21