Buscando una función que se ajuste a una curva tipo sigmoide

Question

Estoy buscando una función para ajustar curvas de tipo sigmoide, a partir de puntos de datos experimentales.

El modelo (la función) no importa, no tiene que ser físicamente relevante, sólo quiero poder calcular y a partir de cualquier x. Y no quiero extrapolar entre dos puntos.

He aquí un ejemplo:

Y aquí están los datos brutos correspondientes:

| X             | Y              |
|---------------|----------------|
| 0             | 0              |
| 1,6366666667  | -12,2012787905 |
| 3,2733333333  | -13,7833876716 |
| 4,91          | -10,5943208589 |
| 6,5466666667  | -1,3584575518  |
| 8,1833333333  | 8,1590423167   |
| 9,82          | 13,8827937482  |
| 10,4746666667 | 18,4965880076  |
| 11,4566666667 | 42,1205206106  |
| 11,784        | 45,0528073182  |
| 12,4386666667 | 76,8150755186  |
| 13,0933333333 | 80,0883540997  |
| 14,73         | 89,7784173678  |
| 16,3666666667 | 98,8113459392  |
| 19,64         | 104,104366506  |
| 22,9133333333 | 105,9929585305 |
| 26,1866666667 | 94,0070414695  |

¿Tienes alguna idea? Mi problema es que los datos van por debajo de 0 en algunos puntos.

EDITAR:

A algunos les molesta el último punto. Para aclarar: al final de la curva, debería haber una meseta. El último punto es sólo un poco de error. Probablemente lo eliminaré de los datos cuando empiece a ajustarlo.

Answer 1

Creo que las splines de suavizado con pequeños grados de libertad servirían. Aquí hay un ejemplo en R:

El código R:

txt <- "| 0             | 0              |
| 1.6366666667  | -12.2012787905 |
| 3.2733333333  | -13.7833876716 |
| 4.91          | -10.5943208589 |
| 6.5466666667  | -1.3584575518  |
| 8.1833333333  | 8.1590423167   |
| 9.82          | 13.8827937482  |
| 10.4746666667 | 18.4965880076  |
| 11.4566666667 | 42.1205206106  |
| 11.784        | 45.0528073182  |
| 12.4386666667 | 76.8150755186  |
| 13.0933333333 | 80.0883540997  |
| 14.73         | 89.7784173678  |
| 16.3666666667 | 98.8113459392  |
| 19.64         | 104.104366506  |
| 22.9133333333 | 105.9929585305 |
| 26.1866666667 | 94.0070414695  |"

dat <- read.table(text=txt, sep="|")[,2:3]
names(dat) <- c("x", "y")
plot(dat$y~dat$x, pch = 19, xlab = "x", ylab = "y", main = "Smoothing Splines with Varying df")

spl3 <- smooth.spline(x = dat$x, y = dat$y, df = 3)
lines(spl3, col = 2)

spl8 <- smooth.spline(x = dat$x, y = dat$y, df = 8)
lines(spl8, col = 4)

legend("topleft", c("df = 3", "df = 8"), col = c(2,4), bty = "n", lty = 1)

Answer 2

Para ajustar una función de tipo sigmoide de forma no paramétrica, podríamos utilizar un spline monótona . Esto se implementa en el paquete R (todos los paquetes R aquí referenciados están en CRAN) splines2 . Tomaré prestado algún código R de la respuesta de @Chaconne, y lo modificaré para mis necesidades.

splines2 ofrece las funciones mSpline implementando M-splines, que es una base de splines no negativa en todas partes (en el intervalo donde se define), y iSpline la integral de la base M-spline. Las últimas son entonces monótonas crecientes, por lo que podemos ajustar una función creciente utilizándolas como base spline de regresión, y ajustar un modelo lineal con restricciones en los coeficientes para que sean no negativos. Esto último se implementa de forma sencilla mediante el paquete R colf "optimización restringida en funciones lineales". Los ajustes parecen:

El código R utilizado:

library(splines2) # includes monotone splines,  M-splines,  I-splines.
library(colf) # constrained optimization on linear functions

 txt <- "| 0             | 0              |
    | 1.6366666667  | -12.2012787905 |
    | 3.2733333333  | -13.7833876716 |
    | 4.91          | -10.5943208589 |
    | 6.5466666667  | -1.3584575518  |
    | 8.1833333333  | 8.1590423167   |
    | 9.82          | 13.8827937482  |
    | 10.4746666667 | 18.4965880076  |
    | 11.4566666667 | 42.1205206106  |
    | 11.784        | 45.0528073182  |
    | 12.4386666667 | 76.8150755186  |
    | 13.0933333333 | 80.0883540997  |
    | 14.73         | 89.7784173678  |
    | 16.3666666667 | 98.8113459392  |
    | 19.64         | 104.104366506  |
    | 22.9133333333 | 105.9929585305 |
    | 26.1866666667 | 94.0070414695  |"

    dat <- read.table(text=txt, sep="|")[,2:3]
names(dat) <- c("x", "y")
plot(dat$y ~ dat$x, pch = 19, xlab = "x", ylab = "y", main = "Monotone Splines with Varying df")

Imod_df_4  <-  colf_nls(y ~ 1 + iSpline(x, df=4), data=dat, lower=c(-Inf, rep(0, 4)), control=nls.control(maxiter=1000, tol=1e-09, minFactor=1/2048) )
lines(dat$x, fitted(Imod_df_4), col="blue")

Imod_df_6  <-  colf_nls(y ~ 1 + iSpline(x, df=6), data=dat, lower=c(-Inf, rep(0, 6)), control=nls.control(maxiter=1000, tol=1e-09, minFactor=1/2048) )
lines(dat$x, fitted(Imod_df_6), col="orange")

Imod_df_8  <-  colf_nls(y ~ 1 + iSpline(x, df=8), data=dat, lower=c(-Inf, rep(0, 8)), control=nls.control(maxiter=1000, tol=1e-09, minFactor=1/2048) )
lines(dat$x, fitted(Imod_df_8), col="red") 

EDIT

Las restricciones monótonas en una spline son un caso especial de splines de forma restringida y ahora hay uno (de hecho varios) paquetes de R que los implementan simplificando su uso. Haré el ejemplo anterior de nuevo, con uno de esos paquetes. El código R está abajo, usando los datos como se leyó arriba:

library(cgam)

mod_cgam0 <- cgam(y ~ 1+s.incr(x), data=dat, family=gaussian)
summary(mod_cgam0)
Call:
cgam(formula = y ~ 1 + s.incr(x), family = gaussian, data = dat)

Coefficients:
            Estimate  StdErr t.value   p.value    
(Intercept)  43.4925  2.7748  15.674 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 102.2557)

Null deviance:  33749.25  on 16  degrees of freedom 
Residual deviance:  1636.091  on 12.5  observed degrees of freedom 

Approximate significance of smooth terms: 
          edf mixture.of.Beta   p.value    
s.incr(x)   3          0.9515 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
CIC:  7.6873 

De este modo, los nudos (y los grados de libertad) se han seleccionado automáticamente. Para fijar el número de grados de libertad utilice:

mod_cgam1 <- cgam(y ~ 1+s.incr(x, numknots=5), data=dat, family=gaussian)

Un documento de presentación de cgam es aquí (arxiv) .

Answer 3

La curva que muestras parece más bien una función cúbica, $ax^3+bx^2+cx+d$ , ya que los extremos giran hacia arriba y hacia abajo, en lugar de extenderse de forma plana/horizontal. O algo así, hecho con una línea de tendencia polinómica en Excel:

Pero si quiere que los extremos se extiendan horizontalmente, hay muchas distribuciones de probabilidad CDF sigmoidales para elegir. Las preguntas que hay que hacerse para elegir la distribución más adecuada son:

1. ¿Cuál es el mecanismo/razón subyacente de una curva de forma sigmoidal?
2. ¿Qué flexibilidad de forma debe tener? ¿Cuántos grados de libertad? Esto dependerá del número de puntos de datos, ya que se quiere evitar el sobreajuste. Pero también, ¿qué características varían y qué características permanecen constantes? ¿La media? ¿La varianza (dispersión)? ¿La asimetría (el lado inclinado)? ¿Kurtosis (colas)?
3. Entonces puedes buscar la forma adecuada en esta lista de la wikipedia ( https://en.wikipedia.org/wiki/List_of_probability_distributions ), o afine su pregunta con más detalles para obtener la mejor respuesta.
4. También hay distribuciones de 4 y 5 parámetros basadas en la función logit con mucha más flexibilidad en la forma, pero de nuevo, debe evitarlas a menos que tenga muchos puntos de datos.

Y PS. Deberías nunca añadir o eliminar selectivamente puntos de datos para el ajuste - ¡BAD BOY!

Answer 4

Puede utilizar la función sigmoid() del paquete {pracma} en R. La función ajustará una curva sigmoidal a un vector numérico. Si no te importa qué función se ajusta a los datos, yo recomendaría la función gam() del paquete {mgcv} en R. Ésta ajusta una función de suavizado a los datos utilizando la regresión spline (el valor por defecto es la placa fina, pero puedes consultar la documentación para otros tipos). Utilizando gam(), como con cualquier ajuste de modelo no paramétrico, no podrá predecir y a partir de x con ninguna fiabilidad fuera del rango de x en su conjunto de datos, ya que las predicciones simplemente seguirán la dirección de la "última pendiente" de la curva, pero por su pregunta parece que no le preocupa eso. Espero que esto ayude.

Answer 5

Junto con las otras sugerencias, una curva de crecimiento de Gompertz también se ajustaría a estos datos. Esto es lo que dice Wikipedia al respecto:

Una curva de Gompertz o función de Gompertz, llamada así por Benjamin Gompertz es una función sigmoidea. Es un tipo de modelo matemático para una serie temporal de tiempo, en el que el crecimiento es más lento al principio y al final de un período de tiempo. ( https://en.wikipedia.org/wiki/Gompertz_function )

La clave es la función sigmoidea. Esta es la fórmula para $y$ en función de $t$ :

$$y(t)=a \exp[-b\exp(-ct)],$$

donde $a$ es una asíntota, ya que $\lim_{t \to \infty} a \exp[-b \exp(-ct)]= a \exp(0)=a$ ;

$b>0$ establece el desplazamiento a lo largo del $x$ -Eje (traslada el gráfico a la izquierda o a la derecha)

$c > 0$ establece la tasa de crecimiento ( $y$ de la escala).

Aquí $\exp(1) = e$ es el número de Euler $2.71828...$ .