Este team de trabajo ha estado mucho tiempo investigando respuestas a tus preguntas, te compartimos la resolución de modo que esperamos resultarte de gran ayuda.
Solución:
El resumen ejecutivo
Puede utilizar el incorporado Ellipsoid funcionan directamente con la media y la covarianza calculadas. Para un 95% de confianza, use:
Ellipsoid[mean, cov Quantile[ChiSquareDistribution[2], 0.95]]
Esa expresión devuelve un Ellipsoid objeto que puede visualizar como un Epilog a un ListPlot, o como argumento para Graphics (más formato a continuación).
Los elipsoides para otros valores críticos comunes se pueden obtener de la misma manera. Tenga en cuenta los diferentes multiplicadores para cov:
90% : Ellipsoid[mean, cov Quantile[ChiSquareDistribution[2], 0.9]]
95% : Ellipsoid[mean, cov Quantile[ChiSquareDistribution[2], 0.95]]
99% : Ellipsoid[mean, cov Quantile[ChiSquareDistribution[2], 0.99]]
Una respuesta más detallada
Primero permítanme comenzar mencionando que una matriz de covarianza debe ser simétrica. Te falta un dígito en cov[[1,2]] que hace que tu original no sea simétrico; Supongo que es un error tipográfico y usaré la versión simétrica a continuación:
mean = 0.968479, 0.020717;
cov = 0.0000797131, 0.000069929, 0.0000699293, 0.00174702
La forma más fácil de generar un elipsoide con la ubicación y alineación correctas dada su distribución es alimentar la media y la covarianza directamente al Ellipsoid función, simplemente como Ellipsoid[mean, cov]. La resultante Ellipsoid es una primitiva gráfica, por lo que puede trazarse sobre sus datos usando, por ejemplo, Epilog o Graphics.
Para obtener un ejemplo práctico, permítanos generar y graficar algunos puntos aleatorios de su distribución, asumiendo que es normal.
SeedRandom[1];
sampledata = RandomVariate[MultinormalDistribution[mean, cov], 2500];
ListPlot[
sampledata,
PlotRange -> All, PlotRangePadding -> Scaled[0.05],
AspectRatio -> 1, Axes -> None, Frame -> True,
Epilog -> Opacity[0], EdgeForm[Thick, Red], Ellipsoid[mean, cov]
]
Como puede ver, sin embargo, un Ellipsoid es decir, “una covarianza de ancho”, como el que trazamos, contiene solo una pequeña fracción de los puntos muestreados (solo aproximadamente el 40% de los puntos, ver más abajo). En su lugar, solicitó un elipsoide que contenga el 95% de los puntos de su distribución. Necesitamos una mayorEllipsoidpara eso: pero ¿qué tan ancho?
Podemos averiguarlo usando la función de distribución de probabilidad (PDF) para su distribución multivariante. Podemos integrar el PDF sobre una región elipsoidal paramétrica para calcular qué fracción de las muestras cae dentro de esa región. Consideremos un bidimensional MultinormalDistribution, con promedio cero y covarianza expresada por sigmax^2, rho sigmax sigmay, rho sigmax sigmay, sigmay^2, dónde sigmax y sigmay son las desviaciones estándar asociadas con cada una de las dos variables independientes, y rho es el coeficiente de correlación entre las dos variables. Las desviaciones estándar son números positivos y 0 <= rho <= 1.
Aquí calculamos una expresión para la fracción de puntos que se encuentran dentro de una elipse bidimensional centrada alrededor de cero (la media de esta distribución) y “norte amplias covarianzas “(observe las n factor en el Ellipsoiddescriptor de. La integración se lleva a cabo sobre la región definida por ese “norte-elipse ancha “.
gencovar = sigmax^2, rho sigmax sigmay, rho sigmax sigmay, sigmay^2;
Assuming[
n > 0, sigmax > 0, sigmay > 0, 0 <= rho < 1,
Simplify[
Integrate[
PDF[MultinormalDistribution[0, 0, gencovar], x, y],
x, y [Element] Ellipsoid[0, 0, n gencovar]
]
]
]
(* 1-E^(-n/2) *)
Ahora tabulemos el valor de esa expresión para unos pocos n.
TableForm[#, TableAlignments -> Right, Top] &@
Table[[email protected]
Esto significa que nuestro elipsoide de "ancho único" original contenía sólo el 39% de las muestras; para lograr una inclusión del 95%, necesitamos 6x anchoEllipsoid.
Tracemos eso para su distribución original (observe el factor 6x muy importante en el Ellipsoid definición):
ListPlot[
sampledata,
PlotRange -> All, PlotRangePadding -> Scaled[0.05],
AspectRatio -> 1, Axes -> None, Frame -> True,
Epilog -> Opacity[0], EdgeForm[Thick, Red], Ellipsoid[mean, 6 cov]
]
Finalmente, también podemos confirmar esto contando explícitamente las muestras que caen dentro de este Ellipsoid. La expresion Select[sampledata, RegionMember[Ellipsoid[mean, n cov]]] nos permite seleccionar esas muestras en sampledata que se encuentran dentro de la región geométrica definida por "6x de ancho" Ellipsoid[mean, 6 cov].
[email protected]@Select[sampledata, RegionMember[Ellipsoid[mean, 6 cov]]] / [email protected]
(* 0.9544 *)
Como se esperaba de nuestros cálculos anteriores, aproximadamente el 95% de los puntos residen dentro del elipsoide 6x definimos.
Para mayor claridad, [email protected][sampledata, RegionMember[Ellipsoid[mean, n cov]]] es el número de puntos que se encuentran dentro del Ellipsoid región. [email protected] es el número total de puntos en nuestra muestra. N ¿Existe para obtener una respuesta numérica aproximada, en lugar de una simbólica?
Por qué no usar EllipsoidQuantile?
EllipsoidQuantile es una función disponible en el MultivariateStatistics paquete que fue mencionado por @Michael E2 en un comentario como una posible solución a este problema. EllipsoidQuantile[dataset, q] devuelve un elipsoide centrado en la media de dataset y escalado para contener una fracción q de El dataset.
A simple vista, esto parecería exactamente lo que pidió el OP, pero esta función se comporta de una manera sutilmente diferente. En mi entendimiento, trata dataset como el La población entera, en lugar de como una muestra de una población más grande. Si la muestra disponible es grande y representa bien a la población, los resultados de los dos métodos serán esencialmente indistinguibles. Sin embargo, los dos métodos darán diferentes resultados para muestras más pequeñas y altos niveles de confianzaq.
También tengo dos razones más prácticas para no usar la función EllipsoidQuantile. Primero, es difícil aplicar estilos a su salida, aunque nunca he podido precisar por qué exactamente. Además, algunas definiciones en MultivariateStatistics parecen sombrear las definiciones de las funciones que desde entonces han pasado al estado integrado, p. ej. PrincipalComponents y MultinormalDistribution, así que prefiero no cargar el paquete a menos que sea absolutamente necesario.
Aquí hay algunos Manipulate código que permite comparar los dos resultados en el sampledata Genere arriba, y un ejemplo de cuando los dos enfoques difieren considerablemente.
Needs["MultivariateStatistics`"]
Manipulate[
Show[
ListPlot[
sampledata[[1 ;; ;; every]],
Axes -> None, Frame -> True,
Epilog ->
Inset[Style["alpha = " <> [email protected], FontSize -> 14], Scaled[0.85, 0.9]]
],
Graphics[
(* using EllipsoidQuantile *)
EllipsoidQuantile[sampledata[[1 ;; ;; every]], alpha],
(* Using the Ellipsoid method outlined above *)
Opacity[0], EdgeForm[Thick, Red],
Ellipsoid[
[email protected][[1 ;; ;; every]],
Evaluate[n /. [email protected]@Solve[1 - E^(-n/2) == alpha, n]] [email protected][[1 ;; ;; every]]
]
]
, PlotRange -> 0.93, 1.01, -0.17, 0.24
],
(* Manipulate variables *)
alpha, 0.95, "[Alpha] value", 0.9, 0.99, 0.01, Appearance -> "Open",
every, 100, "points", 1 -> "All", 10 -> "250", 20 -> "125", 50 -> "50", 100 -> "25", 250 -> "10", ControlType -> SetterBar
]
Una salida de muestra destaca la diferencia: en rojo está el Ellipsoid resultado, y en negro el EllipsoidQuantile producción:
Código de trabajo con Ellipsoid y Eigensystem
<< MultivariateStatistics`
mean = 0, 0;
cov = 1, 0.7, 0.7, 1;
sampledata = RandomVariate[MultinormalDistribution[mean, cov], 2500];
Show[
ListPlot[sampledata, PlotRange -> All],
Graphics[
Ellipsoid[0, 0, Sqrt[2*6*Eigensystem[cov][[1]]],
Eigensystem[cov][[2]]]],
PlotRange -> -5, 5, -5, 5, AspectRatio -> 1
]
F
Nos encantaría que puedieras dar recomendación a este ensayo si te valió la pena.