Newton-Raphson on the Cauchy Likelihood
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 1000
library(shiny)
### ---------- Cauchy location model ------------------------------------------
## Negative log-likelihood of a location theta for data x (the constant
## n*log(pi) is dropped), its gradient g = L' and its curvature J = L''.
nll <- function(th, x)
vapply(th, function(t) sum(log1p((x - t)^2)), numeric(1))
gfun <- function(th, x)
vapply(th, function(t) -sum(2 * (x - t) / (1 + (x - t)^2)), numeric(1))
jfun <- function(th, x)
vapply(th, function(t) sum(2 * (1 - (x - t)^2) / (1 + (x - t)^2)^2),
numeric(1))
### ---------- robustification of the curvature -------------------------------
## The Newton step is -g / J. Where J <= 0 (the tails) it points uphill, and
## where J is close to 0 it is enormous. Two repairs replace J by J~ > 0:
## ridge : J~ = J + tau, with tau >= 0 the smallest value giving
## J~ >= rho * n/2 (Levenberg-Marquardt; in one dimension this
## is the same as flooring the curvature at rho * n/2)
## scoring : J~ = n/2, the expected information of n Cauchy observations
adjust_curv <- function(J, method, n, rho) {
switch(method,
none = J,
ridge = J + max(0, rho * n / 2 - J),
scoring = n / 2)
}
### ---------- Newton-Raphson path --------------------------------------------
## Returns the iterates with L, g, J, the curvature actually used (J~), the
## proposed step and the fraction s of it that was accepted. Stops on
## convergence, or when the iteration breaks down.
newton_path <- function(th0, x, maxit = 25, halve = FALSE, robust = "none",
rho = 0.2, tol = 1e-9) {
n <- length(x)
cn <- c("theta", "nll", "g", "J", "Jt", "step", "s")
P <- matrix(NA_real_, maxit + 1, length(cn), dimnames = list(NULL, cn))
fill <- function(i, th) {
J <- jfun(th, x)
P[i, 1:5] <<- c(th, nll(th, x), gfun(th, x), J,
adjust_curv(J, robust, n, rho))
}
th <- th0; fill(1, th)
n_ok <- 1; note <- "maximum number of iterations reached"
for (i in seq_len(maxit)) {
g <- P[i, 3]; Jt <- P[i, 5]
if (abs(g) < 1e-12) { note <- "converged: the gradient is zero"; break }
if (!is.finite(Jt) || Jt == 0) {
note <- "the curvature is zero: the Newton step is undefined"
break
}
step <- -g / Jt
P[i, 6] <- step
s <- 1
## step halving: shrink s until the Armijo condition
## L(theta + s*step) <= L(theta) + c1 * s * g * step holds
if (halve) {
if (g * step >= 0) {
note <- paste("the step points uphill (the curvature is not",
"positive), so halving cannot help: stopped")
break
}
f0 <- P[i, 2]
while (s > 1e-12) {
fn <- nll(th + s * step, x)
if (is.finite(fn) && fn <= f0 + 1e-4 * s * g * step) break
s <- s / 2
}
if (s <= 1e-12) { note <- "halving found no decrease: stopped"; break }
}
th <- th + s * step
if (!is.finite(th)) { note <- "the iterate is no longer finite: diverged"; break }
P[i, 7] <- s
fill(i + 1, th)
n_ok <- i + 1
if (abs(th) > 1e6) { note <- "diverged: the iterate has run off to infinity"; break }
if (abs(s * step) < tol) { note <- "converged"; break }
}
list(P = P[seq_len(n_ok), , drop = FALSE], n = n_ok, note = note)
}
### ---------- UI --------------------------------------------------------------
## an icon-only playback button; the title shows as a tooltip
ctrl_btn <- function(id, icon_name, title, class = "btn-default")
actionButton(id, NULL, icon = icon(icon_name), title = title, class = class,
style = "flex: 1; padding: 6px 0;")
ui <- fluidPage(
titlePanel("Shinylive App for Newton-Raphson on the Cauchy Likelihood"),
sidebarLayout(
sidebarPanel(
width = 4,
radioButtons("data", "data",
c("one observation, x = 0" = "single",
"random Cauchy sample" = "sample")),
conditionalPanel(
"input.data == 'sample'",
div(style = "display: flex; gap: 10px;",
div(style = "flex: 1;",
numericInput("n", "sample size n", value = 5, min = 1,
max = 50, step = 1, width = "100%")),
div(style = "flex: 1;",
numericInput("seed", "random seed", value = 1,
step = 1, width = "100%")))),
numericInput("th0", HTML("starting value θ<sub>0</sub>"),
value = 1.2, step = 0.1, width = "100%"),
helpText("Or click inside either plot to choose a new starting ",
"value."),
selectInput("halve", "step length",
c("full Newton step" = "none",
"step halving (Armijo)" = "halve")),
selectInput("robust", "robustification of the Hessian",
c("none (observed curvature J)" = "none",
"ridge: J + τ, so that J~ ≥ ε" = "ridge",
"Fisher scoring: J~ = n/2" = "scoring")),
conditionalPanel(
"input.robust == 'ridge'",
sliderInput("rho", HTML("floor ε as a fraction of n/2"),
min = 0.02, max = 1, value = 0.2, step = 0.02)),
sliderInput("W", "half-width of the plotted window", min = 2,
max = 30, value = 6, step = 1),
checkboxInput("compare", "overlay the plain Newton-Raphson path",
TRUE),
checkboxInput("autoplay", "autoplay when the start changes", TRUE),
sliderInput("iter", "iteration", min = 0, max = 1, value = 0,
step = 1),
div(style = "display: flex; gap: 4px;",
ctrl_btn("first", "backward-fast", "Back to the start"),
ctrl_btn("back", "backward-step", "One step back"),
ctrl_btn("play", "play", "Play", class = "btn-primary"),
ctrl_btn("pause", "pause", "Pause"),
ctrl_btn("step", "forward-step", "One step forward"),
ctrl_btn("toend", "forward-fast", "To the end")),
hr(),
helpText("With one observation at 0, plain Newton-Raphson ",
"converges only from starts with |θ₀| < ",
"1/√3 ≈ 0.58. Try 0.5, then 0.8 and 1.2, ",
"first with both safeguards off, then with each in turn: ",
"halving alone stalls once the curvature turns ",
"negative (|θ| > 1), because the step then points ",
"uphill; a positive curvature J~ repairs the direction, ",
"and halving then tames its length.")
),
mainPanel(
width = 8,
plotOutput("pL", click = "clickL", height = "260px"),
plotOutput("pg", click = "clickg", height = "260px"),
hr(),
div(style = "max-width: 640px;",
tableOutput("tab"),
htmlOutput("note"))
)
)
)
### ---------- server ----------------------------------------------------------
server <- function(input, output, session) {
playing <- reactiveVal(FALSE)
play_ms <- 900
## clicking either plot moves the starting value
set_start <- function(cl) {
if (!is.null(cl) && is.finite(cl$x))
updateNumericInput(session, "th0", value = round(cl$x, 2))
}
observeEvent(input$clickL, set_start(input$clickL))
observeEvent(input$clickg, set_start(input$clickg))
x <- reactive({
if (input$data == "single") return(0)
req(input$n, input$seed)
set.seed(input$seed)
rcauchy(max(1, round(input$n)), location = 0)
})
## the plotted window is centred on the median of the data
win <- reactive(median(x()) + c(-1, 1) * input$W)
start <- reactive({ req(is.finite(input$th0)); input$th0 })
modified <- reactive(input$halve != "none" || input$robust != "none")
path <- reactive(newton_path(start(), x(), halve = input$halve == "halve",
robust = input$robust,
rho = if (is.null(input$rho)) 0.2 else input$rho))
## the plain (unsafeguarded) path, shown for comparison
ref <- reactive({
if (!isTRUE(input$compare) || !modified()) return(NULL)
newton_path(start(), x())
})
## curves on a fixed grid over the window
grid <- reactive({
th <- seq(win()[1], win()[2], length.out = 600)
list(th = th, L = nll(th, x()), g = gfun(th, x()))
})
## current iteration, clamped to the length of the path
k <- reactive(min(as.integer(input$iter), path()$n - 1))
## a new path restarts the display at iteration 0
observeEvent(path(), {
updateSliderInput(session, "iter", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
go_to <- function(i) {
updateSliderInput(session, "iter",
value = max(0, min(i, path()$n - 1)))
}
observeEvent(input$play, {
if (k() >= path()$n - 1) go_to(0)
if (path()$n > 1) playing(TRUE)
})
observeEvent(input$pause, playing(FALSE))
observeEvent(input$step, { playing(FALSE); go_to(k() + 1) })
observeEvent(input$back, { playing(FALSE); go_to(k() - 1) })
observeEvent(input$first, { playing(FALSE); go_to(0) })
observeEvent(input$toend, { playing(FALSE); go_to(path()$n - 1) })
observe({
if (!isTRUE(playing())) return()
kk <- isolate(k())
if (kk >= isolate(path()$n) - 1) { playing(FALSE); return() }
invalidateLater(play_ms, session)
isolate(go_to(kk + 1))
})
in_win <- function(v) v >= win()[1] & v <= win()[2]
## the marks shared by both plots: data, the grid minimum, the plain path
base_marks <- function() {
rug(x(), col = "darkgreen", lwd = 2, ticksize = 0.05)
}
output$pL <- renderPlot({
G <- grid(); P <- path()$P; kk <- k() + 1L; RP <- ref()
par(mar = c(4, 4.5, 2, 1))
plot(G$th, G$L, type = "l", lwd = 2, col = "grey30",
xlab = expression(theta), ylab = expression(L(theta)),
main = "negative log-likelihood and the Newton iterates")
base_marks()
i0 <- which.min(G$L)
points(G$th[i0], G$L[i0], pch = 4, col = "darkgreen", lwd = 2, cex = 1.3)
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
Q <- Q[in_win(Q[, 1]), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "grey45", lty = 2, lwd = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "grey85", col = "grey45", cex = 0.9)
}
Q <- P[1:kk, , drop = FALSE]; Q <- Q[in_win(Q[, 1]), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "steelblue", lwd = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
if (in_win(cur[1])) {
points(cur[1], cur[2], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (theta = %.3g) is outside the plotted window",
cur[1]), side = 3, line = -1.5, col = "firebrick",
cex = 0.9)
}
leg <- c("current iterate", "minimum in the window", "data")
legend("top", bg = "#ffffffcc", box.col = NA, cex = 0.85, ncol = 3,
pch = c(21, 4, NA), pt.bg = c("firebrick", NA, NA),
lwd = c(NA, 2, 2), lty = c(NA, NA, 1),
col = c("black", "darkgreen", "darkgreen"), legend = leg)
})
output$pg <- renderPlot({
G <- grid(); P <- path()$P; kk <- k() + 1L; RP <- ref()
par(mar = c(4, 4.5, 2, 1))
plot(G$th, G$g, type = "l", lwd = 2, col = "grey30",
xlab = expression(theta), ylab = expression(g(theta)),
main = "gradient and the tangent (Newton) steps")
abline(h = 0, lty = 3)
base_marks()
## each step follows the line through (theta_k, g_k) whose slope is
## the curvature that was used, J~, down to the horizontal axis
draw_steps <- function(Q, col, lwd, lty) {
m <- nrow(Q)
if (m < 2) return()
for (i in seq_len(m - 1)) {
segments(Q[i, 1], Q[i, 3], Q[i, 1] + Q[i, 6], 0,
col = col, lwd = lwd, lty = lty)
if (isTRUE(Q[i, 7] < 1)) # halved: mark where the full step lands
points(Q[i, 1] + Q[i, 6], 0, pch = 4, col = col, cex = 0.9)
else
segments(Q[i + 1, 1], 0, Q[i + 1, 1], Q[i + 1, 3],
col = "grey40", lty = 3)
}
}
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
draw_steps(Q, "grey65", 1.5, 2)
Q <- Q[in_win(Q[, 1]), , drop = FALSE]
points(Q[, 1], Q[, 3], pch = 21, bg = "grey85", col = "grey45", cex = 0.9)
}
draw_steps(P[1:kk, , drop = FALSE], "steelblue", 2, 1)
Q <- P[1:kk, , drop = FALSE]; Q <- Q[in_win(Q[, 1]), , drop = FALSE]
points(Q[, 1], Q[, 3], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
if (in_win(cur[1])) {
points(cur[1], cur[3], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (theta = %.3g) is outside the plotted window",
cur[1]), side = 3, line = -1.5, col = "firebrick",
cex = 0.9)
}
})
## significant digits throughout; a diverging run reaches values that would
## otherwise be far too wide
fs <- function(v) ifelse(!is.finite(v), "", formatC(v, digits = 4, format = "g"))
## why the run ended, shown as one line under the table
output$note <- renderUI({
R <- path(); RP <- ref()
if (k() + 1 < R$n) return(NULL)
HTML(paste0("<i>", R$note,
if (!is.null(RP)) paste0("; plain Newton-Raphson: ", RP$note),
"</i>"))
})
output$tab <- renderTable({
P <- path()$P[seq_len(k() + 1), , drop = FALSE]
d <- data.frame(iteration = as.character(seq_len(nrow(P)) - 1),
theta = fs(P[, 1]), `L(theta)` = fs(P[, 2]),
g = fs(P[, 3]), J = fs(P[, 4]), `J~ used` = fs(P[, 5]),
`step kept` = ifelse(is.na(P[, 7]), "", fs(P[, 7])),
check.names = FALSE, stringsAsFactors = FALSE)
tail(d, 6)
}, rownames = FALSE, align = "r", width = "100%")
}
shinyApp(ui, server)About the app
Shows how the Newton-Raphson step diverges on the univariate Cauchy location likelihood, and how step halving and a robustified Hessian each repair it. The upper plot is the negative log-likelihood \(L(\theta) = \sum_i \log\{1 + (x_i - \theta)^2\}\) with the iterates marked on it; the lower plot is its gradient \(g = L'\), where each Newton step is drawn as the line through \((\theta_k, g_k)\) whose slope is the curvature used, followed down to the horizontal axis.
Click in either plot (or type a value) to choose the starting value \(\theta_0\), then use the playback buttons to step through the iterations or play the path as an animation. With one observation at \(0\) the gradient \(g(\theta) = 2\theta/(1+\theta^2)\) flattens in the tails and the curvature \(J = L''\) turns negative for \(|\theta| > 1\), so the tangent points away from the minimum and the iterates run off to infinity; plain Newton-Raphson converges only from starts with \(|\theta_0| < 1/\sqrt{3}\). Choose a random sample instead to see the same behaviour on a likelihood with several local minima.
The two safeguards act on different parts of the step \(-g/\tilde J\). Step halving shortens a step until the Armijo condition holds, so it tames the length but cannot fix a step that points uphill. Robustifying the Hessian replaces \(J\) by a positive \(\tilde J\), which fixes the direction: ridge adds the smallest \(\tau \ge 0\) that makes \(J + \tau \ge \varepsilon\), and Fisher scoring uses the expected information \(n/2\). In one dimension ridging and flooring the curvature coincide. Switching on either keeps the plain iteration visible as a dashed grey path for comparison; the details of each iteration, and the reason the run ended, are reported under the plots.
NoteR source for this app
library(shiny)
### ---------- Cauchy location model ------------------------------------------
## Negative log-likelihood of a location theta for data x (the constant
## n*log(pi) is dropped), its gradient g = L' and its curvature J = L''.
nll <- function(th, x)
vapply(th, function(t) sum(log1p((x - t)^2)), numeric(1))
gfun <- function(th, x)
vapply(th, function(t) -sum(2 * (x - t) / (1 + (x - t)^2)), numeric(1))
jfun <- function(th, x)
vapply(th, function(t) sum(2 * (1 - (x - t)^2) / (1 + (x - t)^2)^2),
numeric(1))
### ---------- robustification of the curvature -------------------------------
## The Newton step is -g / J. Where J <= 0 (the tails) it points uphill, and
## where J is close to 0 it is enormous. Two repairs replace J by J~ > 0:
## ridge : J~ = J + tau, with tau >= 0 the smallest value giving
## J~ >= rho * n/2 (Levenberg-Marquardt; in one dimension this
## is the same as flooring the curvature at rho * n/2)
## scoring : J~ = n/2, the expected information of n Cauchy observations
adjust_curv <- function(J, method, n, rho) {
switch(method,
none = J,
ridge = J + max(0, rho * n / 2 - J),
scoring = n / 2)
}
### ---------- Newton-Raphson path --------------------------------------------
## Returns the iterates with L, g, J, the curvature actually used (J~), the
## proposed step and the fraction s of it that was accepted. Stops on
## convergence, or when the iteration breaks down.
newton_path <- function(th0, x, maxit = 25, halve = FALSE, robust = "none",
rho = 0.2, tol = 1e-9) {
n <- length(x)
cn <- c("theta", "nll", "g", "J", "Jt", "step", "s")
P <- matrix(NA_real_, maxit + 1, length(cn), dimnames = list(NULL, cn))
fill <- function(i, th) {
J <- jfun(th, x)
P[i, 1:5] <<- c(th, nll(th, x), gfun(th, x), J,
adjust_curv(J, robust, n, rho))
}
th <- th0; fill(1, th)
n_ok <- 1; note <- "maximum number of iterations reached"
for (i in seq_len(maxit)) {
g <- P[i, 3]; Jt <- P[i, 5]
if (abs(g) < 1e-12) { note <- "converged: the gradient is zero"; break }
if (!is.finite(Jt) || Jt == 0) {
note <- "the curvature is zero: the Newton step is undefined"
break
}
step <- -g / Jt
P[i, 6] <- step
s <- 1
## step halving: shrink s until the Armijo condition
## L(theta + s*step) <= L(theta) + c1 * s * g * step holds
if (halve) {
if (g * step >= 0) {
note <- paste("the step points uphill (the curvature is not",
"positive), so halving cannot help: stopped")
break
}
f0 <- P[i, 2]
while (s > 1e-12) {
fn <- nll(th + s * step, x)
if (is.finite(fn) && fn <= f0 + 1e-4 * s * g * step) break
s <- s / 2
}
if (s <= 1e-12) { note <- "halving found no decrease: stopped"; break }
}
th <- th + s * step
if (!is.finite(th)) { note <- "the iterate is no longer finite: diverged"; break }
P[i, 7] <- s
fill(i + 1, th)
n_ok <- i + 1
if (abs(th) > 1e6) { note <- "diverged: the iterate has run off to infinity"; break }
if (abs(s * step) < tol) { note <- "converged"; break }
}
list(P = P[seq_len(n_ok), , drop = FALSE], n = n_ok, note = note)
}
### ---------- UI --------------------------------------------------------------
## an icon-only playback button; the title shows as a tooltip
ctrl_btn <- function(id, icon_name, title, class = "btn-default")
actionButton(id, NULL, icon = icon(icon_name), title = title, class = class,
style = "flex: 1; padding: 6px 0;")
ui <- fluidPage(
titlePanel("Shinylive App for Newton-Raphson on the Cauchy Likelihood"),
sidebarLayout(
sidebarPanel(
width = 4,
radioButtons("data", "data",
c("one observation, x = 0" = "single",
"random Cauchy sample" = "sample")),
conditionalPanel(
"input.data == 'sample'",
div(style = "display: flex; gap: 10px;",
div(style = "flex: 1;",
numericInput("n", "sample size n", value = 5, min = 1,
max = 50, step = 1, width = "100%")),
div(style = "flex: 1;",
numericInput("seed", "random seed", value = 1,
step = 1, width = "100%")))),
numericInput("th0", HTML("starting value θ<sub>0</sub>"),
value = 1.2, step = 0.1, width = "100%"),
helpText("Or click inside either plot to choose a new starting ",
"value."),
selectInput("halve", "step length",
c("full Newton step" = "none",
"step halving (Armijo)" = "halve")),
selectInput("robust", "robustification of the Hessian",
c("none (observed curvature J)" = "none",
"ridge: J + τ, so that J~ ≥ ε" = "ridge",
"Fisher scoring: J~ = n/2" = "scoring")),
conditionalPanel(
"input.robust == 'ridge'",
sliderInput("rho", HTML("floor ε as a fraction of n/2"),
min = 0.02, max = 1, value = 0.2, step = 0.02)),
sliderInput("W", "half-width of the plotted window", min = 2,
max = 30, value = 6, step = 1),
checkboxInput("compare", "overlay the plain Newton-Raphson path",
TRUE),
checkboxInput("autoplay", "autoplay when the start changes", TRUE),
sliderInput("iter", "iteration", min = 0, max = 1, value = 0,
step = 1),
div(style = "display: flex; gap: 4px;",
ctrl_btn("first", "backward-fast", "Back to the start"),
ctrl_btn("back", "backward-step", "One step back"),
ctrl_btn("play", "play", "Play", class = "btn-primary"),
ctrl_btn("pause", "pause", "Pause"),
ctrl_btn("step", "forward-step", "One step forward"),
ctrl_btn("toend", "forward-fast", "To the end")),
hr(),
helpText("With one observation at 0, plain Newton-Raphson ",
"converges only from starts with |θ₀| < ",
"1/√3 ≈ 0.58. Try 0.5, then 0.8 and 1.2, ",
"first with both safeguards off, then with each in turn: ",
"halving alone stalls once the curvature turns ",
"negative (|θ| > 1), because the step then points ",
"uphill; a positive curvature J~ repairs the direction, ",
"and halving then tames its length.")
),
mainPanel(
width = 8,
plotOutput("pL", click = "clickL", height = "260px"),
plotOutput("pg", click = "clickg", height = "260px"),
hr(),
div(style = "max-width: 640px;",
tableOutput("tab"),
htmlOutput("note"))
)
)
)
### ---------- server ----------------------------------------------------------
server <- function(input, output, session) {
playing <- reactiveVal(FALSE)
play_ms <- 900
## clicking either plot moves the starting value
set_start <- function(cl) {
if (!is.null(cl) && is.finite(cl$x))
updateNumericInput(session, "th0", value = round(cl$x, 2))
}
observeEvent(input$clickL, set_start(input$clickL))
observeEvent(input$clickg, set_start(input$clickg))
x <- reactive({
if (input$data == "single") return(0)
req(input$n, input$seed)
set.seed(input$seed)
rcauchy(max(1, round(input$n)), location = 0)
})
## the plotted window is centred on the median of the data
win <- reactive(median(x()) + c(-1, 1) * input$W)
start <- reactive({ req(is.finite(input$th0)); input$th0 })
modified <- reactive(input$halve != "none" || input$robust != "none")
path <- reactive(newton_path(start(), x(), halve = input$halve == "halve",
robust = input$robust,
rho = if (is.null(input$rho)) 0.2 else input$rho))
## the plain (unsafeguarded) path, shown for comparison
ref <- reactive({
if (!isTRUE(input$compare) || !modified()) return(NULL)
newton_path(start(), x())
})
## curves on a fixed grid over the window
grid <- reactive({
th <- seq(win()[1], win()[2], length.out = 600)
list(th = th, L = nll(th, x()), g = gfun(th, x()))
})
## current iteration, clamped to the length of the path
k <- reactive(min(as.integer(input$iter), path()$n - 1))
## a new path restarts the display at iteration 0
observeEvent(path(), {
updateSliderInput(session, "iter", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
go_to <- function(i) {
updateSliderInput(session, "iter",
value = max(0, min(i, path()$n - 1)))
}
observeEvent(input$play, {
if (k() >= path()$n - 1) go_to(0)
if (path()$n > 1) playing(TRUE)
})
observeEvent(input$pause, playing(FALSE))
observeEvent(input$step, { playing(FALSE); go_to(k() + 1) })
observeEvent(input$back, { playing(FALSE); go_to(k() - 1) })
observeEvent(input$first, { playing(FALSE); go_to(0) })
observeEvent(input$toend, { playing(FALSE); go_to(path()$n - 1) })
observe({
if (!isTRUE(playing())) return()
kk <- isolate(k())
if (kk >= isolate(path()$n) - 1) { playing(FALSE); return() }
invalidateLater(play_ms, session)
isolate(go_to(kk + 1))
})
in_win <- function(v) v >= win()[1] & v <= win()[2]
## the marks shared by both plots: data, the grid minimum, the plain path
base_marks <- function() {
rug(x(), col = "darkgreen", lwd = 2, ticksize = 0.05)
}
output$pL <- renderPlot({
G <- grid(); P <- path()$P; kk <- k() + 1L; RP <- ref()
par(mar = c(4, 4.5, 2, 1))
plot(G$th, G$L, type = "l", lwd = 2, col = "grey30",
xlab = expression(theta), ylab = expression(L(theta)),
main = "negative log-likelihood and the Newton iterates")
base_marks()
i0 <- which.min(G$L)
points(G$th[i0], G$L[i0], pch = 4, col = "darkgreen", lwd = 2, cex = 1.3)
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
Q <- Q[in_win(Q[, 1]), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "grey45", lty = 2, lwd = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "grey85", col = "grey45", cex = 0.9)
}
Q <- P[1:kk, , drop = FALSE]; Q <- Q[in_win(Q[, 1]), , drop = FALSE]
if (nrow(Q) > 1) lines(Q[, 1], Q[, 2], col = "steelblue", lwd = 2)
points(Q[, 1], Q[, 2], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
if (in_win(cur[1])) {
points(cur[1], cur[2], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (theta = %.3g) is outside the plotted window",
cur[1]), side = 3, line = -1.5, col = "firebrick",
cex = 0.9)
}
leg <- c("current iterate", "minimum in the window", "data")
legend("top", bg = "#ffffffcc", box.col = NA, cex = 0.85, ncol = 3,
pch = c(21, 4, NA), pt.bg = c("firebrick", NA, NA),
lwd = c(NA, 2, 2), lty = c(NA, NA, 1),
col = c("black", "darkgreen", "darkgreen"), legend = leg)
})
output$pg <- renderPlot({
G <- grid(); P <- path()$P; kk <- k() + 1L; RP <- ref()
par(mar = c(4, 4.5, 2, 1))
plot(G$th, G$g, type = "l", lwd = 2, col = "grey30",
xlab = expression(theta), ylab = expression(g(theta)),
main = "gradient and the tangent (Newton) steps")
abline(h = 0, lty = 3)
base_marks()
## each step follows the line through (theta_k, g_k) whose slope is
## the curvature that was used, J~, down to the horizontal axis
draw_steps <- function(Q, col, lwd, lty) {
m <- nrow(Q)
if (m < 2) return()
for (i in seq_len(m - 1)) {
segments(Q[i, 1], Q[i, 3], Q[i, 1] + Q[i, 6], 0,
col = col, lwd = lwd, lty = lty)
if (isTRUE(Q[i, 7] < 1)) # halved: mark where the full step lands
points(Q[i, 1] + Q[i, 6], 0, pch = 4, col = col, cex = 0.9)
else
segments(Q[i + 1, 1], 0, Q[i + 1, 1], Q[i + 1, 3],
col = "grey40", lty = 3)
}
}
if (!is.null(RP)) {
Q <- RP$P[seq_len(min(kk, RP$n)), , drop = FALSE]
draw_steps(Q, "grey65", 1.5, 2)
Q <- Q[in_win(Q[, 1]), , drop = FALSE]
points(Q[, 1], Q[, 3], pch = 21, bg = "grey85", col = "grey45", cex = 0.9)
}
draw_steps(P[1:kk, , drop = FALSE], "steelblue", 2, 1)
Q <- P[1:kk, , drop = FALSE]; Q <- Q[in_win(Q[, 1]), , drop = FALSE]
points(Q[, 1], Q[, 3], pch = 21, bg = "steelblue", cex = 1.1)
cur <- P[kk, ]
if (in_win(cur[1])) {
points(cur[1], cur[3], pch = 21, bg = "firebrick", cex = 1.8)
} else {
mtext(sprintf("current iterate (theta = %.3g) is outside the plotted window",
cur[1]), side = 3, line = -1.5, col = "firebrick",
cex = 0.9)
}
})
## significant digits throughout; a diverging run reaches values that would
## otherwise be far too wide
fs <- function(v) ifelse(!is.finite(v), "", formatC(v, digits = 4, format = "g"))
## why the run ended, shown as one line under the table
output$note <- renderUI({
R <- path(); RP <- ref()
if (k() + 1 < R$n) return(NULL)
HTML(paste0("<i>", R$note,
if (!is.null(RP)) paste0("; plain Newton-Raphson: ", RP$note),
"</i>"))
})
output$tab <- renderTable({
P <- path()$P[seq_len(k() + 1), , drop = FALSE]
d <- data.frame(iteration = as.character(seq_len(nrow(P)) - 1),
theta = fs(P[, 1]), `L(theta)` = fs(P[, 2]),
g = fs(P[, 3]), J = fs(P[, 4]), `J~ used` = fs(P[, 5]),
`step kept` = ifelse(is.na(P[, 7]), "", fs(P[, 7])),
check.names = FALSE, stringsAsFactors = FALSE)
tail(d, 6)
}, rownames = FALSE, align = "r", width = "100%")
}
shinyApp(ui, server)This app accompanies Optimization for Maximum Likelihood Estimation in the book.