# Loading datasets
data(simple)
head(simple)
################################################
# Time delay estimation via profile likelihood #
################################################
###### The entire profile likelihood values on the grid of values of the time delay.
# Cubic microlensing model
ti1 <- simple[, 1]
ti2 <- simple[, 1]^2
ti3 <- simple[, 1]^3
ss <- lm(simple[, 4] - simple[, 2]~ ti1 + ti2 + ti3)
initial <- c(mean(simple[, 2]), log(0.01), log(200), ss$coefficients)
delta.uniform.range <- c(0, 100)
grid <- seq(0, 100, by = 0.1)
# grid interval "by = 0.1" is recommended.
### Running the following codes takes more time than CRAN policy
### Please type the following lines without "#" to run the function and to see the results
# logprof <- entirelogprofilelikelihood(data = simple, grid = grid,
# initial = initial, data.flux = FALSE,
# delta.uniform.range = delta.uniform.range, micro = 3)
# plot(grid, logprof, type = "l",
# xlab = expression(bold(Delta)),
# ylab = expression(bold(paste("log L"[prof], "(", Delta, ")"))))
# prof <- exp(logprof - max(logprof)) # normalization
# plot(grid, prof, type = "l",
# xlab = expression(bold(Delta)),
# ylab = expression(bold(paste("L"[prof], "(", Delta, ")"))))
Run the code above in your browser using DataLab