Hi I'm currently doing coding to simulate data using inverse method. Im using parallel exponential model where I let the lambda=b0+b1x. My simulation is based on survival analysis.
#generate data
gen <- function(n,lambda,b0,b1){
set.seed(1)
u <- runif(n,0,1)
c1 <- rexp(n,lambda)
x <- rnorm(n,0,1)
t1 = -log(1 - sqrt(u) ) / (b0 + b1*x) #inverse method
c <- 1*(t1 < c1)
t = pmin(t1, c1)
data1 <- data.frame(x, t, t1, c1, c)
return(data1)
}
data2 <- gen(20,0.01,2,4)
data2
x = data2$x
t = data2$t
xsum = sum(x)
tsum = sum(t)
The problem is that when run the second coding below, it won't show my mle for b0 and b1
#Likelihood
library(maxLik)
LLF <- function(para){
set.seed(1)
b0 = para[1]
b1 = para[2]
n = 1
z1 = (n*log(2)) + (n*log(b0+b1*xsum)) - ((b0+b1*xsum)*tsum) + (n*log(1-exp((-(b0 + b1*xsum)*tsum))))
return(z1)
}
mle <- maxLik(LLF, start = c(2,4))
The problem is you assigned n=1 in the LLF. Since we usually maximize the parameters given the entire data, n should be equal to number of observations. If you update this info, your mle will converge. For example,
n<-nrow(data2)
#Likelihood
library(maxLik)
LLF <- function(para){
set.seed(1)
b0 = para[1]
b1 = para[2]
#n = 1
z1 = (n*log(2)) + (n*log(b0+b1*xsum)) - ((b0+b1*xsum)*tsum) + (n*log(1-exp((-(b0 + b1*xsum)*tsum))))
return(z1)
}
mle <- maxLik(LLF, start = c(2,4))
summary(mle)
Maximum Likelihood estimation
Newton-Raphson maximisation, 3 iterations
Return code 1: gradient close to zero
Log-Likelihood: -22.7055
2 free parameters
Estimates:
Estimate Std. error t value Pr(> t)
[1,] 1.986 NA NA NA
[2,] 3.986 NA NA NA