Skip to content

Levenberg–Marquardt en R

De Wikipedia:

In mathematics and computing, the Levenberg–Marquardt algorithm (LMA or just LM), also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. These minimization problems arise especially in least squares curve fitting.

The LMA is used in many software applications for solving generic curve-fitting problems. However, as with many fitting algorithms, the LMA finds only a local minimum, which is not necessarily the global minimum. The LMA interpolates between the Gauss–Newton algorithm (GNA) and the method of gradient descent.

Hay una cosa del Graphpad medio misteriosa que tiene implementación en R, pero faltaba cómo usarlo para ajustar "el histograma de unos datos a una normal".

Adapté el códido de este post en R-Bloggers: https://www.r-bloggers.com/a-better-nls/

Bueno:

20200421-170917.png

En negro está la PDF estimada de los datos con density() y en rojo el ajuste de Levenberg–Marquardt.

Código

library(minpack.lm)
library(tidyverse)

#### Simular datos ####

mean.1 <- 10
mean.2 <- 15

d <- c(rnorm(1000, mean = mean.1, sd = 1),
  rnorm(100, mean = mean.2, sd = 1)
  )

#### Graficar datos simulados ####

ggplot(data.frame(d = d)) + geom_density(aes(x = d))

#### Estimar PDF y graficarla ####

# Deberia dar lo mismo que el geom_density de antes
d.d <- density(d)
d.d.df <- data.frame(x = d.d$x, density=d.d$y)

ggplot(dd.d.df) + geom_point(aes(x = x, y = density))

#### Fittear la PDF con Levenberg–Marquardt ####

d.LM <- nlsLM(density ~ ( 1/(Sigma*sqrt(2*pi)) ) * exp(1)^(-0.5*((x-Mu)/Sigma)^2), 
                   data = d.d.df, 
                   start = list(Mu = 9, Sigma = 1)
)

# d.LM <- nlsLM(density ~ dnorm(x = x, sd = 1, mean = Mu), 
#                    data = d.d.df, 
#                    start = list(Mu = 9)
# )

#### Extraer coeficientes y graficar la PDF fitteada sobre la original ####

d.LM.coef <- coef(d.LM)

# Agregar fitteo al dataframe
d.d.df$density.fit <- dnorm(d.d.df$x, mean = d.LM.coef["Mu"], sd = d.LM.coef["Sigma"])

ggplot(d.d.df) + 
  geom_point(aes(x = x, y = density)) + 
  geom_line(aes(x = x, y = density.fit), color = "red")