Sample points from along a ridge This "dents" the likelihood surface by reflecting points better than a threshold back across the threshold (think of taking a hollow plastic model of a mountain and punching the top so it's a volcano). It then uses essentially a Metropolis-Hastings walk to wander around the new rim. It adjusts the proposal width so that it samples points around the desired likelihood. This is better than using the curvature at the maximum likelihood estimate since it can actually sample points in case the assumptions of the curvature method do not hold. It is better than varying one parameter at a time while holding others constant because that could miss ridges: if I am fitting 5=x+y, and get a point estimate of (3,2), the reality is that there are an infinite range of values of x and y that will sum to 5, but if I hold x constant it looks like y is estimated very precisely. Of course, one could just fully embrace the Metropolis-Hastings lifestyle and use a full Bayesian approach.

Source: R/dentist.R

dent_walk.Rd

While running, it will display the current range of likelihoods in the desired range (by default, the best negative log likelihood + 2 negative log likelihood units) and the parameter values falling in that range. If things are working well, the range of values will stabilize during a search.

dent_walk(
  par,
  fn,
  best_neglnL,
  delta = 2,
  nsteps = 1000,
  print_freq = 50,
  lower_bound = 0,
  upper_bound = Inf,
  adjust_width_interval = 100,
  badval = 1e+09,
  sd_vector = NULL,
  debug = FALSE,
  restart_after = 50,
  ...
)

Arguments

par: Starting parameter vector, generally at the optimum. If named, the vector names are used to label output parameters.
fn: The likelihood function, assumed to return negative log likelihoods
best_neglnL: The negative log likelihood at the optimum; other values will be greater than this.
delta: How far from the optimal negative log likelihood to focus samples
nsteps: How many steps to take in the analysis
print_freq: Output progress every print_freq steps.
lower_bound: Minimum parameter values to try. One for all or a vector of the length of par.
upper_bound: Maximum parameter values to try. One for all or a vector of the length of par.
adjust_width_interval: When to try automatically adjusting proposal widths
badval: Bad negative log likelihood to return if a non-finite likelihood is returned
sd_vector: Vector of the standard deviations to use for proposals. Generated automatically if NULL
debug: If TRUE, prints out much more information during a run
restart_after: Sometimes the search can get stuck outside the good region but still accept moves. After this many steps without being inside the good region, restart from one of the past good points
...: Other arguments to fn.

Value

A dentist object containing results, the data.frame of negative log likelihoods and the parameters associated with them; acceptances, the vector of whether a proposed move was accepted each step; best_neglnL, the best value passed into the analysis; delta, the desired offset; all_ranges, a summary of the results.

Details

The algorithm tunes: if it is moving too far away from the desired likelihoods, it will decrease the proposal width; if it staying in areas better than the desired likelihood, it will increase the proposal width. It will also expand the proposal width for parameters where the extreme values still appear good enough to try to find out the full range for these values.

In general, the idea of this is not to give you a pleasingly narrow range of possible values – it is to try to find the actual uncertainty, including finding any ridges that would not be seen in univariate space.

Examples

# Univariate case
sims <- stats::rnorm(100, mean=17)
possible_means <- seq(from=16, to=18, length.out=100) # for optimize

# Make sure we have a function that takes in a parameters vector, other arguments if needed,
# and returns the negative log likelihood
dnorm_to_run <- function(par, sims) {
  return(-sum(stats::dnorm(x=sims, mean=par, log=TRUE)))
}

optimized_results <- stats::optimize(dnorm_to_run,interval=range(possible_means), 
  sims=sims, maximum=FALSE)
best_par <- optimized_results$minimum
names(best_par) <- "mean"
best_neglnL <- optimized_results$objective

dented_results <- dent_walk(par=best_par, fn=dnorm_to_run, best_neglnL=best_neglnL, sims=sims)
#> [1] "Done replicate 50"
#> [1] "CI of values (the 6 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 17.07602
#> [2,] 144.4501 17.29485
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 100"
#> [1] "CI of values (the 10 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.95104
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 150"
#> [1] "CI of values (the 18 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 200"
#> [1] "CI of values (the 24 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 250"
#> [1] "CI of values (the 32 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 300"
#> [1] "CI of values (the 42 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 350"
#> [1] "CI of values (the 48 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 400"
#> [1] "CI of values (the 55 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 450"
#> [1] "CI of values (the 58 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 500"
#> [1] "CI of values (the 68 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 550"
#> [1] "CI of values (the 77 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 600"
#> [1] "CI of values (the 90 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 650"
#> [1] "CI of values (the 99 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 700"
#> [1] "CI of values (the 106 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 750"
#> [1] "CI of values (the 116 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 800"
#> [1] "CI of values (the 128 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 850"
#> [1] "CI of values (the 142 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 900"
#> [1] "CI of values (the 154 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "Done replicate 950"
#> [1] "CI of values (the 164 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.91344
#> [2,] 144.4501 17.29485
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 1000"
#> [1] "CI of values (the 174 replicates within 2 neglnL of the optimum)"
#>        neglnL     mean
#> [1,] 142.6182 16.90843
#> [2,] 144.5195 17.29485
plot(dented_results)


# Multivariate case
sims <- stats::rlnorm(100, meanlog=1, sdlog=3)

dlnorm_to_run <- function(par, sims) {
  return(-sum(stats::dlnorm(sims, meanlog=par[1], sdlog=par[2], log=TRUE)))
}

optimized_results <- stats::optim(c(meanlog=.9, sdlog=2.9), dlnorm_to_run, sims=sims)
best_par <- optimized_results$par
best_neglnL <- optimized_results$value

dented_results <- dent_walk(par=best_par, fn=dlnorm_to_run, best_neglnL=best_neglnL, sims=sims)
#> [1] "Done replicate 50"
#> [1] "CI of values (the 20 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.9657083 2.853762
#> [2,] 355.3208 1.3999136 3.610972
#> [1] "increasing proposal width for all parameters"
#> [1] "Done replicate 100"
#> [1] "CI of values (the 29 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.9657083 2.853762
#> [2,] 355.3450 1.5449791 3.610972
#> [1] "Done replicate 150"
#> [1] "CI of values (the 45 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.8222937 2.853762
#> [2,] 355.3450 1.5449791 3.610972
#> [1] "increasing proposal width for some parameters"
#> [1] "Done replicate 200"
#> [1] "CI of values (the 58 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.8222937 2.853762
#> [2,] 355.3701 1.5539521 3.610972
#> [1] "Done replicate 250"
#> [1] "CI of values (the 68 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.8222937 2.853762
#> [2,] 355.3701 1.5539521 3.610972
#> [1] "Done replicate 300"
#> [1] "CI of values (the 82 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.8222937 2.853762
#> [2,] 355.3701 1.5539521 3.610972
#> [1] "Done replicate 350"
#> [1] "CI of values (the 95 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.6897976 2.853762
#> [2,] 355.3701 1.5539521 3.610972
#> [1] "increasing proposal width for all parameters"
#> [1] "Done replicate 400"
#> [1] "CI of values (the 105 replicates within 2 neglnL of the optimum)"
#>        neglnL  meanlog    sdlog
#> [1,] 353.3707 0.638462 2.784415
#> [2,] 355.3701 1.553952 3.648617
#> [1] "Done replicate 450"
#> [1] "CI of values (the 114 replicates within 2 neglnL of the optimum)"
#>        neglnL  meanlog    sdlog
#> [1,] 353.3707 0.395771 2.784415
#> [2,] 355.3701 1.553952 3.648617
#> [1] "Done replicate 500"
#> [1] "CI of values (the 123 replicates within 2 neglnL of the optimum)"
#>        neglnL  meanlog    sdlog
#> [1,] 353.3707 0.395771 2.784415
#> [2,] 355.3701 1.553952 3.648617
#> [1] "Done replicate 550"
#> [1] "CI of values (the 131 replicates within 2 neglnL of the optimum)"
#>        neglnL  meanlog    sdlog
#> [1,] 353.3707 0.395771 2.784415
#> [2,] 355.3701 1.553952 3.648617
#> [1] "Done replicate 600"
#> [1] "CI of values (the 139 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 650"
#> [1] "CI of values (the 146 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "decreasing proposal width for all parameters"
#> [1] "Done replicate 700"
#> [1] "CI of values (the 153 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 750"
#> [1] "CI of values (the 165 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 800"
#> [1] "CI of values (the 170 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 850"
#> [1] "CI of values (the 183 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.784415
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 900"
#> [1] "CI of values (the 189 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.757891
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 950"
#> [1] "CI of values (the 199 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.757891
#> [2,] 355.3701 1.5539521 3.648617
#> [1] "Done replicate 1000"
#> [1] "CI of values (the 211 replicates within 2 neglnL of the optimum)"
#>        neglnL   meanlog    sdlog
#> [1,] 353.3707 0.3408285 2.757891
#> [2,] 355.3701 1.5539521 3.648617
plot(dented_results)