4 votos

N-sigma curvas para un no-lineal de mínimos cuadrados el ajuste de curva

Estoy usando python's scipy.optimizar.curve_fit rutina (que utiliza una no-lineal de mínimos cuadrados) para adaptarse a una función exponencial de la forma:

f(x) = a * exp(b*x) + c

para un conjunto de datos. El resultado se parece a esto:

enter image description here

donde los triángulos negros son el conjunto de datos y la curva azul es la f(x) ajustado por la rutina. La salida del proceso incluye el valor óptimo para cada parámetro a, b, c junto con una matriz de covarianza de donde puedo obtener la desviación estándar para cada uno de ellos, es decir, obtengo:

$a = p_a \pm e_a$

$b = p_b \pm e_b$

$c = p_c \pm e_c$

donde $p_i$ son los valores óptimos y $e_i$ sus desviaciones estándar.

Lo que quiero es generar el N-$\sigma$ curvas (donde N es un entero) para el ajuste. Como ejemplo, he aquí un artículo donde los autores muestran un similar exponencial de mejor ajuste (aunque una descomposición uno), junto con sus $1\sigma$ superior e inferior de las curvas (véase la descripción):

enter image description here

Este sitio proporciona un ejemplo de cómo hacer esto, pero no estoy seguro de que el procedimiento es completamente correcto.

¿Cómo puedo obtener correctamente esas curvas con los datos que tengo disponible?

5voto

AdamSane Puntos 1825

Uno de los grandes problemas que veo con mínimos cuadrados ordinarios es la varianza está claro que no es constante, sino que aumenta a medida que $x$ no; yo estaría inclinado a considerar en un modelo de propagación es la relativa a la media (posiblemente proporcional a ella), lo que sugiere una generalizada no lineales del modelo o tal vez el uso de una transformación (como un registro-transform) para estabilizar la varianza.

Sin tener que lidiar con ese problema, su varianza-covarianza y el error estándar de las estimaciones sesgadas.

Puesto que usted quiere tener estimaciones de las funciones de tales cantidades, abordar este problema va a ser central para la obtención de un utilizable respuesta.


En respuesta a los comentarios:

Los errores estándar de las estimaciones de los parámetros y el residual de la desviación estándar se basa en el supuesto de varianza constante. Este supuesto es violado claramente. Considerar si (a ojo) me dibujar barras verticales destinadas a cubrir alrededor del 99.algo por ciento de los locales de difusión en cada uno de los diferentes valores de x:

enter image description here

Si los cálculos están basados en que la igualdad de la asunción (como la costumbre de los errores estándar son), entonces es probable que sea de poco valor.

Los autores en el papel que su segunda parcela no viene (una página de referencia hubiera sido bueno por cierto - yo no soy un astrónomo, así que tiene poca idea de lo que está pasando en el resto de él) parecen haber utilizado un ajuste donde propagación no es constante para el y escala de su parcela, y de hecho parece aumentar de una manera consistente con sus barras de error, por lo que no está claro que ellos tienen necesariamente un problema.

El blog de enlace para no ver como la trama del papel; parece ser el uso de diferentes supuestos. Cosas que se ven superficialmente similar no puede ser.

Si la propagación se ve más o menos constante en la escala logarítmica (y que no hay real 0), entonces usted podría encajar un no lineal de mínimos cuadrados del modelo de la forma $E(y) = \log(a e^{bx}+c)$, y luego convertirlo a un ajuste en la unstransformed escala, pero tenga en cuenta que las expectativas en las dos escalas no corresponden directamente. Si usted tiene cerca de la normalidad en el registro de la escala que usted todavía puede obtener un estimado de la expectativa en la escala original. En su defecto puede que se quede con una determinada distribución de la asunción, o una aproximación de Taylor a la expectativa.


Varios tipos de intervalos

Existen varios tipos de intervalos de personas posible que desee colocar en una parcela. Aquí puedo hablar de tres:

$\;$ 1. un pointwise (confianza) en el intervalo (condicional error estándar de la media). Este es un intervalo donde la curva (la media condicional), para cada valor de x. Es decir, especificar un determinado $x$, y tomar un corte vertical de allí. Podemos producir una estimación del error estándar de la media, en que $x$; en la repetición de una serie de dos tal error estándar de los intervalos de* (a partir de diferentes datos) en un determinado $x$ sería de esperar para incluir a la población en ese $x$ sobre el 95% del tiempo si ciertas suposiciones estaban satisfechos.

Esto parece ser algo así como lo que se intentó en el post del blog que usted menciona en su último enlace, pero, fundamentalmente, mirando el código de Python parecen estar ignorando el efecto del parámetro de covarianzas, lo que haría que el resultado incorrecto.

$\;$ *(No decimos "$2\sigma$" a pesar de que, desde el $\sigma$ convencionalmente ser una cantidad de la población, y este es un ejemplo de estimación de algo)

$\;$ 2. Un total de intervalo para la media. Esto sería todavía un intervalo para la media, sino una basada en la cobertura global de las probabilidades. (No creo que usted está después de esto.)

$\;$ 3. Un (pointwise) intervalo de predicción. Este sería un intervalo con una probabilidad especificada de la inclusión de una nueva observación en cada $x$. Se incluye tanto la incertidumbre de los parámetros (tanto de varianzas y covarianzas de las estimaciones de los parámetros) y la observación de "ruido" (basado en una muestra de la estimación de lo que yo llamaría $\sigma$). Esto es lo que a primera vista parece estar sucediendo en la segunda trama, pero no estoy 100% seguro.

Usted debe aclarar qué tipo de intervalo que usted necesita, y lo que exactamente quieres decir con "$\sigma$".


Puesto que no estás muy familiarizado con la costumbre de los errores estándar y los intervalos en los modelos lineales, vamos a cubrir los. (El cálculo del intervalo de abajo son todos pointwise, no global.)

Los errores estándar y los intervalos en los modelos lineales

A. error Estándar de la media, el intervalo de confianza para la media

Para un modelo lineal, $y=X\beta+\varepsilon$, $\text{Var}(\varepsilon)=\sigma^2I$, tenemos una matriz de $C$.

$\text{Var}(C\hat{\beta})=\sigma^2C(X'X)^{-1}C'$.

Por lo tanto, si $c$ es una fila del vector de constantes, $\text{se}(c\hat\beta)=\sigma\sqrt{c(X'X)^{-1}c'}$.

Si la estimación de $\hat{\sigma^2}$ $s^2=\frac{1}{n-p}\sum_i e_i^2$ (donde $p$ es el número de factores incluyendo el término de intersección), a continuación, $\frac{c\hat\beta-c\beta}{\text{se}(c\hat\beta)}$ es fundamental la cantidad* con una distribución t con $n-p$ d.f.

Un $1-\alpha$ intervalo de confianza para una media (en $c=x^*$) es, por tanto, $x* \hat\beta \pm t_{1-\alpha/2} s \sqrt{x^* (X'X)^{-1} {x^*}'}$ o, más sucintamente, $\hat{y} \pm t_{1-\alpha/2} s \sqrt{x^* (X'X)^{-1} {x^*}'}$.

* básicamente un importante cantidad es una función de un parámetro o los parámetros y datos cuya distribución no depende del parámetro(s)

B. la variación del Proceso

Esto es sólo la estimación de la desviación estándar del proceso, $\sigma$, generalmente estimado por $s$. Si supiéramos la verdad (condicional) media en $x^*$, uno esperaría que el 68% de las observaciones dentro de$\sigma$, y alrededor del 95% dentro de las $2\sigma$. El hecho de que $s$ se estima que altera un poco (pero yo no trabajo el punto de dar todos los detalles aquí a menos que usted los necesita, ya que esto no suele ser lo que la gente después).

C. error Estándar de la predicción, el intervalo de predicción de una observación

Este incorpora tanto de los dos anteriores componentes. Deje $y^*$ ser una observación no en nuestra muestra.

Entonces la varianza de $y^{*}-\hat{y}^{*}$$\sigma^2+\sigma^2x^*(X'X)^{-1} {x^*}' = \sigma^2(1+x^* (X'X)^{-1} {x^*}')$, y así el error estándar es $\sigma$ $\sqrt{1+x^* (X'X)^{-1} {x^*}'}$

Ampliación de (ligeramente) la noción de una fundamental a la cantidad, de la misma manera que el intervalo de confianza se puede construir una $1-\alpha$ intervalo de $y^{*}$: $\hat{y}^*$ $\pm$ $t_{1-\alpha/2}$ $s$ $\sqrt{1+x^* (X'X)^{-1} {x^*}'}$


Similar (si es aproximado) los intervalos pueden ser fácilmente construidos en modelos no lineales; por ejemplo, a través de una expansión de Taylor. Si eres capaz de identificar cual (si alguna) de los cálculos anteriores que usted necesita, podemos ser capaces de explorar el modo de enfoque.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X