1 Gradiente Descendente

Para materializar o algoritmo de gradiente descendente (GD), vamos considerar a estimação de parâmetros em modelos de regressão não linear. Primeiramente, modelos não lineares são estimados considerando algoritmos baseados em Newton-Raphson que fazem uso da primeira e segunda derivadas do modelo em relação aos parâmetros. O GD só utiliza-se da primeira derivada.

1.1 Simulação do modelo

# Carregando pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)
#-----------------------------------------------------------------------
# Simulando dados do modelo Michaelis Menten.

# Retorna os valores determinados pelo modelo para a média de Y|x.
eta <- function(x, theta) {
    theta[1] * x/(theta[2] + x)
}

# Valores para os parâmetros.
pars <- c(tha = 10,
          thv = 2)

# Gráfico da função média.
curve(eta(x, pars),
      from = 0,
      to = 20)

# Retorna um rúido normal padrão.
epsilon <- function(x, ...) {
    rnorm(n = length(x), mean = 0, ...)
}

# Simulando do modelo.
set.seed(321)
da <- data.frame(x = sort(runif(n = 50, min = 3, max = 20)))
da$y <- eta(x = da$x, theta = pars) + epsilon(x = da$x, sd = 0.1)

plot(y ~ x, data = da)
curve(eta(x, theta = pars),
      add = TRUE,
      col = "orange")

1.2 Inspeção da função objetivo

Vamos avaliar a função objetivo (função perda ou custo) para um grid de valores dos parâmetros.

# Função perda (custo ou objetivo).
loss <- function(theta, y, x) {
    sum((y - eta(x, theta = theta))^2)
}

loss(theta = pars, y = da$y, x = da$x)
## [1] 0.4824309
# Grid de valores dos parâmetros.
grid <- expand.grid(tha = seq(0, 20, length.out = 41),
                    thv = seq(-3, 8, length.out = 41))
head(grid)
##   tha thv
## 1 0.0  -3
## 2 0.5  -3
## 3 1.0  -3
## 4 1.5  -3
## 5 2.0  -3
## 6 2.5  -3
# Avalia a função objetivo em cada ponto do grid.
grid$rss <- mapply(tha = grid$tha,
                   thv = grid$thv,
                   FUN = function(tha, thv) {
                       loss(theta = c(tha, thv),
                            y = da$y,
                            x = da$x)
                   })

# Define escala de cores para visualização.
colr <- colorRampPalette(rev(brewer.pal(11, "Spectral")),
                         space = "rgb")

# Exibe a superfície da função objetivo (com log para mais detalhes).
levelplot(log(rss) ~ tha + thv,
          data = grid,
          col.regions = colr,
          contour = TRUE) +
    layer(panel.abline(v = pars["tha"],
                       h = pars["thv"],
                       lty = 2))

wireframe(log(rss) ~ tha + thv,
          data = grid,
          drape = TRUE,
          col.regions = colr(101),
          panel.3d.wireframe = wzRfun::panel.3d.contour,
          type = c("bottom"))

1.3 O algoritmo

A função objetivo (custo) para o problema em mãos é \[ L(\theta) = \sum_{i = 1}^{n} (y_{i} - f(x_{i}, \theta))^2, \] em que \(f\) é o modelo considerado para a média de \(Y|x\), ou seja, \[ f(x_{i}, \theta) = \frac{\theta_a x_i }{\theta_v + x_i}, \quad \theta = (\theta_a, \theta_v)'. \]

A derivada primeira da função \(L\) com relação a cada um dos seus argumentos fornece o vetor gradiente. Você pode usar o Wolfran para ganhar tempo no computo dessas derivadas. Por exemplo, passe a expressão para o motor do Wolfran “derivative of (y - a * x/(v + x))^2 wrt a”. Você também pode usar o R para fazer essas derivadas de forma analítica ou numérica. Outra opção é usar o wxMaxima.

#-----------------------------------------------------------------------
# Gradiente.

# Derivadas analíticas com o R.
loss_expr <- expression((y - tha * x/(thv + x))^2)
D(expr = loss_expr, name = "tha")
## -(2 * (x/(thv + x) * (y - tha * x/(thv + x))))
D(expr = loss_expr, name = "thv")
## 2 * (tha * x/(thv + x)^2 * (y - tha * x/(thv + x)))
# Função gradiente.
grad <- function(theta, y, x) {
    tha <- theta[1]
    thv <- theta[2]
    part <-
        cbind(partial_tha =
                  -(2 * (x/(thv + x) * (y - tha * x/(thv + x)))),
              partial_thv =
                  2 * (tha * x/(thv + x)^2 * (y - tha * x/(thv + x))))
    # Derivada da soma = soma das derivadas.
    colSums(part)
}

# Avaliando a função objetivo.
loss(theta = c(10, 2), y = da$y, x = da$x)
## [1] 0.4824309
# Avaliando a função gradiente.
grad(theta = c(10, 2),
     y = da$y,
     x = da$x)
## partial_tha partial_thv 
##   0.5346204  -0.7181287
#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente.

# Tolerância: parar quando a mudança for irrelevante ou exceder número
# máximo de iterações.
tol <- 1E-5
niter <- 10000

# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.001

# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)

# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(tha = 18, thv = -1)
repeat {
    delta <- alpha * grad(theta = th[k, ],
                          y = da$y,
                          x = da$x)
    th[k + 1, ] <- th[k, ] - delta
    # print(th[k, ])
    if (k > niter | all(abs(delta) < tol)) break
    k <- k + 1
}

th <- th[1:k, ]

#-----------------------------------------------------------------------
# Traço do algoritmo.

levelplot(rss ~ tha + thv,
          data = grid,
          contour = TRUE) +
    layer(panel.abline(v = pars["tha"],
                       h = pars["thv"],
                       lty = 2)) +
    layer(panel.points(x = th[, 1],
                       y = th[, 2],
                       type = "o",
                       cex = 0.1))

# Número de iterações.
k
## [1] 2412
# Estimativas.
th[k, ]
## [1] 10.035399  2.056582
# A descida.
lth <- apply(th,
             MARGIN = 1,
             FUN = loss,
             y = da$y,
             x = da$x)
plot(log(lth) ~ seq_along(lth), type = "o", cex = 0.3)

#-----------------------------------------------------------------------
# Faça testes.

# Repita as instruções acima com os valores abaixo e veja o algorítmo
# divergir.
th[k, ] <- c(tha = 10, thv = 3)
alpha <- 0.1

#-----------------------------------------------------------------------
# Ajustando o modelo usando um Gauss-Newton.

n0 <- nls(y ~ tha * x/(thv + x),
          data = da,
          start = c(tha = 18, thv = -1),
          trace = TRUE)
## 7019.99 :  18 -1
## 75.90063 :   9.0302143 -0.2983704
## 4.555524 :  9.495021 1.143315
## 0.5184406 :  9.966329 1.943018
## 0.4720646 :  10.032693  2.052832
## 0.4720599 :  10.033171  2.053785
## 0.4720599 :  10.033167  2.053779
summary(n0)
## 
## Formula: y ~ tha * x/(thv + x)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## tha 10.03317    0.04758  210.87   <2e-16 ***
## thv  2.05378    0.05890   34.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09917 on 48 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 5.96e-08

2 Gradiente Descendente Estocástico

Em termos simples, o gradiente descendente estocástico difere apenas no fato de não usar todos os dados mas sim usar uma fração dos dados. Esse subconjunto dos dados é obtido com reamostragem dos dados originais. A reamostragem é o componente estocástico do algorítmo. Com isso, o traço da função objetivo deixa de ser suave, o que impossibilita o uso de critérios de parada baseados em diferenças de valores consecutivos. Então o mais comum é executar o algoritmo até exceder o número máximo de iterações. O número máximo de iterações deve ser escolhido de forma a garantir suficiente proximidade com o ponto de ótimo da função.

#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente Estocástico.

# Para quando exceder número máximo de iterações.
niter <- 30000

# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.01

# Número de elementos amostrados do conjunto de treinamento.
index <- 1:nrow(da)
nele <- 3

# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)

# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(tha = 18, thv = -1)
while (k < niter) {
    db <- da[sample(index, size = nele), ]
    delta <- alpha * grad(theta = th[k, ],
                          y = db$y,
                          x = db$x)
    th[k + 1, ] <- th[k, ] - delta
    # print(th[k, ])
    k <- k + 1
}

#-----------------------------------------------------------------------
# O traço do algorítmo.

th <- th[1:k, ]
levelplot(rss ~ tha + thv,
          data = grid,
          contour = TRUE) +
    layer(panel.abline(v = pars["tha"],
                       h = pars["thv"],
                       lty = 2)) +
    layer(panel.points(x = th[, 1],
                       y = th[, 2],
                       type = "l"))

# Estimativas.
th[k, ]
## [1] 10.03792  2.04956
# A descida.
lth <- apply(th,
             MARGIN = 1,
             FUN = loss,
             y = da$y,
             x = da$x)
plot(log(lth), type = "o", cex = 0.3)

# Vendo apenas as últimas iterações.
plot((tail(lth, n = 500)), type = "o", cex = 0.3)

3 Gradiente Descendente Boosting

#-----------------------------------------------------------------------
# Gerando dados.

curve(10 + 0.05 * x + 2 * sin(x/10), 0, 100)

# Retorna os valores determinados pelo modelo.
eta <- function(x, beta = c(10, 0.05, 2)) {
    beta[1] + beta[2] * x + beta[3] * sin(x/10)
}

# Simulando do modelo.
set.seed(567)
da <- data.frame(x = sort(runif(200, 0, 100)))
da$y <- eta(da$x) + epsilon(da$x, sd = 1)

# Exibindo os dados.
plot(y ~ x, data = da)
curve(eta, add = TRUE, col = "orange", lwd = 3)

3.1 Exemplo com função “patamar”

A função patamar é uma função por partes definidas sobre o domínio de \(x\) dividido em \(p\) subconjuntos justapostos em que para cada parte a função é de valor constante, ou seja \[ f(x) = \begin{cases} y_1, & x_0 < x \leq x_1 \\ y_2, & x_1 < x \leq x_2 \\ \vdots & \vdots \\ y_p, & x_{p-1} < x \leq x_p. \end{cases} \]

Vamos considerar a versão com \(p = 2\). Dessa forma, a função retornará o valor \(y_1\) para valores \(x < \beta\) e \(y_2\) para valores de \(x \geq \beta\), em que \(\beta\) é o ponto de corte do domínio.

#-----------------------------------------------------------------------
# A função h(x) será uma função "patamar".

# Modelo inicial com predito sendo a média das observações.
m0 <- lm(y ~ 1, data = da)

plot(y ~ x, data = da)
lines(da$x, fitted(m0), type = "o", col = 2)

# Função patamar. Retorna ajustado e resíduo.
h <- function(beta, y, x) {
    dico <- ifelse(x < beta, 0, 1)
    fit <- ave(y, dico, FUN = mean)
    list(dico = dico, fit = fit, res = y - fit)
}

# Função custo.
loss <- function(beta, y, x) {
    fit <- h(beta, y, x)$fit
    sum((y - fit)^2)
}

# Verifica o comportamento da função.
# bseq <- seq(0, 100, 1)
# lseq <- sapply(bseq, FUN = loss, y = da$y, x = da$x)
# plot(lseq ~ bseq, type = "o")

# Usando a optim().
# optim(par = c(50),
#       fn = loss,
#       y = da$y,
#       x = da$x)$par

#-----------------------------------------------------------------------
# O algorítmo.

# Cria uma lista para armazenar os artefatos produzidos.
n <- nrow(da)
k <- 12
resul <- replicate(k,
                   data.frame(fit = numeric(n),
                              res = numeric(n)),
                   simplify = FALSE)

# Inicia a lista com resultados do primeiro ajuste.
resul[[1]]$fit <- fitted(m0)
resul[[1]]$res <- residuals(m0)
bval <- c(NA, numeric(k - 1))

# Executa o gradiente boosting.
for (i in 2:k) {
    # Encontra o beta minimizando a função perda.
    beta <- optimize(f = loss,
                     y = resul[[i - 1]]$res,
                     x = da$x,
                     interval = c(0, 100))$minimum
    # Exibe o progresso.
    cat(i, "\t", beta, "\n")
    bval[i] <- beta
    # Extrai valores ajustados e resíduos.
    hfit <- h(beta,
              y = resul[[i - 1]]$res,
              x = da$x)
    # Atualiza os valores ajustados e resíduos.
    resul[[i]]$fit <- resul[[i - 1]]$fit + hfit$fit
    resul[[i]]$res <- da$y - resul[[i]]$fit
}
## 2     61.12248 
## 3     27.60907 
## 4     69.15174 
## 5     60.47349 
## 6     3.845893 
## 7     23.60687 
## 8     69.28974 
## 9     60.863 
## 10    7.733429 
## 11    26.46976 
## 12    47.80609
#-----------------------------------------------------------------------
# Exibe os resultados.

# Lista replicando a mesma fórmula.
form <- replicate(k, y ~ x)
names(form) <- 1:k

xyplot.list(form,
            data = da,
            type = "n",
            as.table = TRUE) +
    layer({
        i <- which.packet()
        fit <- resul[[i]]$fit
        panel.points(x = x,
                     y = y,
                     col = ifelse(y > fit, "blue", "red"))
        panel.lines(da$x,
                    fit,
                    col = "green2",
                    lwd = 3)
        panel.abline(v = bval[i], lty = 2)
    })

# Converte lista em tabela.
names(resul) <- as.character(1:k)
result <- plyr::ldply(resul, .id = "k")
result$k <- factor(result$k, levels = 1:k)
result$x <- da$x

# O comportamento dos resíduos com o acréscimo de funções.
xyplot(res ~ x | k,
       data = result,
       type = c("n", "smooth"),
       span = 0.3,
       as.table = TRUE) +
    layer(panel.points(x = x, y = y,
                       col = ifelse(y > 0, "blue", "red"))) +
    layer(panel.abline(h = 0, lty = 2))

3.2 Exemplo com funções segmentadas lineares

#-----------------------------------------------------------------------
# Ajustar modelos segmentados.

library(segmented)

# Função que ajusta o modelo segmentado com um ponto de quebra. Já está
# com a função de perda quadrática embutida na segmented().
h <- function(y, x) {
    m <- lm(y ~ x)
    s <- segmented(m, seg.Z = ~x)
    fit <- fitted(s)
    brk <- summary.segmented(s)$psi[1, 2]
    list(brk = brk, fit = fit, res = y - fit)
}

# Verifica o que a função retorna.
str(h(y = da$y, x = da$x))
## List of 3
##  $ brk: num 44.6
##  $ fit: Named num [1:200] 12.1 12.1 12 12 12 ...
##   ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ...
##  $ res: Named num [1:200] -3.17 -3.835 -0.175 -1.134 -2.42 ...
##   ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ...
# Cria uma lista para armazenar os artefatos produzidos.
n <- nrow(da)
k <- 6
resul <- replicate(k,
                   data.frame(fit = numeric(n),
                              res = numeric(n)),
                   simplify = FALSE)

# Inicia a lista com resultados do primeiro ajuste.
resul[[1]]$fit <- fitted(m0)
resul[[1]]$res <- residuals(m0)
bval <- c(NA, numeric(k - 1))

# Executa o gradiente boosting.
for (i in 2:k) {
    # Extrai valores ajustados e resíduos.
    hfit <- h(y = resul[[i - 1]]$res,
              x = da$x)
    # Extrai o ponto de quebra.
    bval[i] <- hfit$brk
    # Atualiza os valores ajustados e resíduos.
    resul[[i]]$fit <- resul[[i - 1]]$fit + hfit$fit
    resul[[i]]$res <- da$y - resul[[i]]$fit
}

#-----------------------------------------------------------------------
# Exibe os resultados.

# Lista replicando a mesma fórmula.
form <- replicate(k, y ~ x)
names(form) <- 1:k

xyplot.list(form,
            data = da,
            type = "n",
            as.table = TRUE) +
    layer({
        i <- which.packet()
        fit <- resul[[i]]$fit
        panel.points(x = x,
                     y = y,
                     col = ifelse(y > fit, "blue", "red"))
        panel.lines(da$x,
                    fit,
                    col = "green2",
                    lwd = 3)
        panel.abline(v = bval[i], lty = 2)
    })

# Converte lista em tabela.
names(resul) <- as.character(1:k)
result <- plyr::ldply(resul, .id = "k")
result$k <- factor(result$k, levels = 1:k)
result$x <- da$x

# O comportamento dos resíduos com o acréscimo de funções.
xyplot(res ~ x | k,
       data = result,
       type = c("n", "smooth"),
       span = 0.3,
       as.table = TRUE) +
    layer(panel.points(x = x, y = y,
                       col = ifelse(y > 0, "blue", "red"))) +
    layer(panel.abline(h = 0, lty = 2))

4 Cuidados com o tunning

#-----------------------------------------------------------------------
# Debugando com problema em uma dimensão.

m0 <- lm(dist ~ speed, data = cars)
coef(m0)
## (Intercept)       speed 
##  -17.579095    3.932409
loss <- function(beta, y, xx) {
    mean((y - (coef(m0)[1] + beta * xx))^2)
}

loss(y = cars$dist, xx = cars$speed, beta = coef(m0)[2])
## [1] 227.0704
deviance(m0)
## [1] 11353.52
D(expression((y - (beta0 + beta * x))^2),
  name = "beta")
## -(2 * (x * (y - (beta0 + beta * x))))
mygrad <- function(y, xx, beta) {
    beta0 <- coef(m0)[1]
    mean(-(2 * (xx * (y - (beta0 + beta * xx)))))
}

library(numDeriv)

numDeriv::grad(func = loss,
               x = coef(m0)[2] + 1,
               y = cars$dist,
               xx = cars$speed)
## [1] 529.12
mygrad(y = cars$dist,
       xx = cars$speed,
       beta = coef(m0)[2] + 1)
## [1] 529.12
# bseq <- seq(0, 8, length.out = 50)
bseq <- seq(-10, 20, length.out = 50)
lseq <- sapply(bseq, FUN = loss, y = cars$dist, x = cars$speed)
plot(lseq ~ bseq, type = "l")

# alpha de 0.001 OK
# alpha de 0.003 OK ZIG-ZAG
# alpha de 0.005 OK

alpha <- 0.0038
n <- 100
bt <- numeric(n)
k <- 1
bt[k] <- -5

while (k < n) {
    gr1 <- numDeriv::grad(func = loss,
                          x = bt[k],
                          y = cars$dist,
                          xx = cars$speed)
    gr2 <- mygrad(y = cars$dist,
                  xx = cars$speed,
                  beta = bt[k])
    l <- loss(y = cars$dist, x = cars$speed, beta = bt[k])
    cat(gr1, "\t", gr2, "\t", bt[k], "\t", l, "\n")
    bt[k + 1] <- bt[k] -
        alpha * gr2
    k <- k + 1
}
## -4726.316     -4726.316   -5      21335.76 
## 4776.68   4776.68     12.96   21788.03 
## -4827.58      -4827.58    -5.191382   22249.98 
## 4879.023      4879.023    13.15342    22721.84 
## -4931.014     -4931.014   -5.386864   23203.8 
## 4983.558      4983.558    13.35099    23696.09 
## -5036.663     -5036.663   -5.586535   24198.93 
## 5090.334      5090.334    13.55279    24712.54 
## -5144.577     -5144.577   -5.790483   25237.15 
## 5199.397      5199.397    13.75891    25773.01 
## -5254.802     -5254.802   -5.998801   26320.34 
## 5310.797      5310.797    13.96945    26879.4 
## -5367.389     -5367.389   -6.211583   27450.45 
## 5424.584      5424.584    14.1845     28033.72 
## -5482.388     -5482.388   -6.428924   28629.49 
## 5540.809      5540.809    14.40415    29238.03 
## -5599.851     -5599.851   -6.650921   29859.61 
## 5659.523      5659.523    14.62851    30494.5 
## -5719.831     -5719.831   -6.877675   31143 
## 5780.782      5780.782    14.85768    31805.39 
## -5842.382     -5842.382   -7.109287   32481.97 
## 5904.638      5904.638    15.09176    33173.05 
## -5967.558     -5967.558   -7.345861   33878.94 
## 6031.148      6031.148    15.33086    34599.95 
## -6095.416     -6095.416   -7.587504   35336.4 
## 6160.369      6160.369    15.57508    36088.64 
## -6226.014     -6226.014   -7.834325   36856.99 
## 6292.358      6292.358    15.82453    37641.81 
## -6359.41      -6359.41    -8.086433   38443.44 
## 6427.176      6427.176    16.07932    39262.25 
## -6495.664     -6495.664   -8.343944   40098.6 
## 6564.881      6564.881    16.33958    40952.87 
## -6634.837     -6634.837   -8.606971   41825.44 
## 6705.538      6705.538    16.60541    42716.71 
## -6776.992     -6776.992   -8.875634   43627.07 
## 6849.207      6849.207    16.87693    44556.94 
## -6922.193     -6922.193   -9.150054   45506.73 
## 6995.955      6995.955    17.15428    46476.87 
## -7070.504     -7070.504   -9.430353   47467.8 
## 7145.848      7145.848    17.43756    48479.96 
## -7221.994     -7221.994   -9.716657   49513.81 
## 7298.951      7298.951    17.72692    50569.8 
## -7376.729     -7376.729   -10.0091    51648.42 
## 7455.335      7455.335    18.02247    52750.15 
## -7534.779     -7534.779   -10.3078    53875.49 
## 7615.07   7615.07     18.32436    55024.94 
## -7696.216     -7696.216   -10.6129    56199.01 
## 7778.227      7778.227    18.63272    57398.24 
## -7861.112     -7861.112   -10.92455   58623.16 
## 7944.88   7944.88     18.94768    59874.33 
## -8029.541     -8029.541   -11.24286   61152.31 
## 8115.103      8115.103    19.26939    62457.66 
## -8201.578     -8201.578   -11.568     63790.99 
## 8288.974      8288.974    19.59799    65152.88 
## -8377.301     -8377.301   -11.90011   66543.95 
## 8466.57   8466.57     19.93364    67964.83 
## -8556.79      -8556.79    -12.23933   69416.15 
## 8647.971      8647.971    20.27647    70898.56 
## -8740.123     -8740.123   -12.58582   72412.74 
## 8833.258      8833.258    20.62665    73959.35 
## -8927.385     -8927.385   -12.93973   75539.11 
## 9022.516      9022.516    20.98434    77152.71 
## -9118.66      -9118.66    -13.30122   78800.88 
## 9215.828      9215.828    21.34968    80484.37 
## -9314.032     -9314.032   -13.67046   82203.93 
## 9413.282      9413.282    21.72286    83960.33 
## -9513.59      -9513.59    -14.04761   85754.36 
## 9614.967      9614.967    22.10403    87586.83 
## -9717.424     -9717.424   -14.43285   89458.56 
## 9820.973      9820.973    22.49336    91370.39 
## -9925.625     -9925.625   -14.82633   93323.19 
## 10031.39      10031.39    22.89104    95317.82 
## -10138.29     -10138.29   -15.22825   97355.19 
## 10246.32      10246.32    23.29724    99436.22 
## -10355.51     -10355.51   -15.63878   101561.8 
## 10465.85      10465.85    23.71214    103733 
## -10577.38     -10577.38   -16.0581    105950.7 
## 10690.09      10690.09    24.13594    108215.8 
## -10804    -10804      -16.48641   110529.6 
## 10919.13      10919.13    24.56881    112892.8 
## -11035.49     -11035.49   -16.92389   115306.8 
## 11153.08      11153.08    25.01095    117772.4 
## -11271.93     -11271.93   -17.37075   120290.9 
## 11392.04      11392.04    25.46257    122863.3 
## -11513.43     -11513.43   -17.82718   125490.9 
## 11636.12      11636.12    25.92387    128174.7 
## -11760.12     -11760.12   -18.29339   130916.1 
## 11885.43      11885.43    26.39505    133716.2 
## -12012.08     -12012.08   -18.76959   136576.2 
## 12140.08      12140.08    26.87632    139497.6 
## -12269.45     -12269.45   -19.256     142481.5 
## 12400.19      12400.19    27.36791    145529.4 
## -12532.33     -12532.33   -19.75282   148642.6 
## 12665.87      12665.87    27.87003    151822.5 
## -12800.84     -12800.84   -20.26029   155070.5 
## 12937.25      12937.25    28.3829     158388.1 
## -13075.11     -13075.11   -20.77863   161776.8 
## 13214.43      13214.43    28.90677    165238.1 
## -13355.25     -13355.25   -21.30808   168773.5
plot(bseq, lseq, type = "l")
lt <- sapply(bt, FUN = loss, y = cars$dist, x = cars$speed)
lines(lt ~ bt, type = "o", col = 2)
abline(v = coef(m0)[2])

head(bt)
## [1] -5.000000 12.960001 -5.191382 13.153422 -5.386864 13.350988
#-----------------------------------------------------------------------
25px