Alicia, miembro de este staff, nos hizo el favor de crear este enunciado ya que domina muy bien este tema.
Solución:
La respuesta de Sjoerd aplica el poder de MathematicaSon herramientas de ajuste de modelos muy generales. Aquí hay una solución de menor tecnología.
Si desea ajustar una distribución gaussiana a un conjunto de datos, puede encontrar su matriz media y de covarianza, y la gaussiana que desea es la que tiene los mismos parámetros. Para los gaussianos, este es en realidad el ajuste óptimo en el sentido de ser el estimador de máxima verosimilitud; para otras distribuciones, esto puede no funcionar igualmente bien, por lo que querrá NonlinearModelFit.
Una arruga es que sus datos no caen a cero sino a algo como Min[data] $ approx 0.0387 $, pero podemos deshacernos de eso fácilmente simplemente restándolo. A continuación, normalizo el array para sumar $ 1 $, de modo que pueda tratarlo como una distribución de probabilidad discreta. (Todo lo que esto realmente hace es permitirme evitar dividir entre Total[data, Infinity] cada vez en lo siguiente.)
min = Min[data];
sum = Total[data - min, Infinity];
p = (data - min)/sum;
Ahora encontramos la media y la covarianza. MathematicaLas funciones de ‘no parecen funcionar con datos ponderados, así que simplemente lanzaremos las nuestras.
m, n = Dimensions[p];
mean = Sum[i, j p[[i, j]], i, 1, m, j, 1, n];
cov = Sum[Outer[Times, i, j - mean, i, j - mean] p[[i, j]], i, 1, m, j, 1, m];
Podemos obtener fácilmente la distribución de probabilidad para el gaussiano con esta media y covarianza. Por supuesto, si queremos hacer coincidir los datos originales, tenemos que reescalarlos y retrocederlos.
f[i_, j_] := PDF[MultinormalDistribution[mean, cov], i, j] // Evaluate;
g[i_, j_] := f[i, j] sum + min;
La coincidencia no es tan mala, aunque los datos tienen dos jorobas donde el gaussiano tiene uno. También puede comparar el ajuste a través de un gráfico de contorno. (Utilice la información sobre herramientas para comparar niveles de contorno).
Show[ListPlot3D[data, PlotStyle -> None],
Plot3D[g[i, j], j, 1, 9, i, 1, 9, MeshStyle -> None,
PlotStyle -> Opacity[0.8]], PlotRange -> All]
Show[ListContourPlot[data, Contours -> Range[0.05, 0.25, 0.025],
ContourShading -> None, ContourStyle -> ColorData[1, 1], InterpolationOrder -> 3],
ContourPlot[g[i, j], j, 1, 9, i, 1, 9, Contours -> Range[0.05, 0.25, 0.025],
ContourShading -> None, ContourStyle -> ColorData[1, 2]]]
Las varianzas a lo largo de los ejes principales son los valores propios de la matriz de covarianza, y las desviaciones estándar (que supongo que llama semiejes) son sus raíces cuadradas.
[email protected][cov]
(* 1.86325, 1.50567 *)
data3D = Flatten[MapIndexed[#2[[1]], #2[[2]], #1 &, data, 2], 1];
fm = NonlinearModelFit[data3D,
a E^(-(((-my + y) Cos[b] - (-mx + x) Sin[b])^2/(2 sy^2)) -
((-mx + x) Cos[b] + (-my + y) Sin[b])^2/(2 sx^2)),
a, 0.1, b, 0, mx, 4.5, my, 4.5, sx, 3, sy, 3,
x, y
]
Show[
ListPlot3D[data3D, PlotStyle -> None, MeshStyle -> Red],
Plot3D[fm["BestFit"], x, 1, 9, y, 1, 9, PlotRange -> All,
PlotStyle -> Opacity[0.9], Mesh -> None]
]
fm["BestFitParameters"]
a -> 0.1871830545, b -> -0.4853901689, mx -> 5.152549499, my -> 5.092511036, sx -> 2.756524919, sy -> 2.072163433
Comenzando con sus datos, trazaría los datos de esta manera:
data3D = Flatten[MapIndexed[#2[[1]], #2[[2]], #1 &, data, 2], 1]
y luego cuéntelo como:
data2 = Flatten[Table[#[[1]], #[[2]], 100 #[[3]]] & /@ data3D, 1];
y luego usa esta función:
FindDistributionParameters[data2,
MultinormalDistribution[a, b, c, d, e, f]]
para encontrar los medios y la cobertura de los datos. Y luego podría hacer otras cosas con la distribución, como crear un objeto de distribución usando:
edist = EstimatedDistribution[data2,
MultinormalDistribution[a, b, c, d, e, f]]
y luego verifique la distribución con los datos, usando:
dtest = DistributionFitTest[data2, edist, "HypothesisTestData"]
y cree una tabla de los resultados de la prueba, usando:
N[dtest["TestDataTable", All]]
o individualmente
AndersonDarlingTest[data2, edist, "TestConclusion"]
CramerVonMisesTest[data2, edist, "TestConclusion"]
JarqueBeraALMTest[data2, "TestConclusion"]
MardiaKurtosisTest[data2, "TestConclusion"]
MardiaSkewnessTest[data2, "TestConclusion"]
PearsonChiSquareTest[data2, edist, "TestConclusion"]
ShapiroWilkTest[data2, "TestConclusion"]
Ahora, no soy un estadístico, y probablemente solo sé suficientes estadísticas para ser peligroso, y no estoy familiarizado con todas las pruebas enumeradas en la tabla, pero sospecho que uno de esos números pequeños es significativo para medir qué tan cerca los datos reales se ajustan a la distribución, es decir, qué tan “gaussiano” es el conjunto de puntos.
También podría trazar los datos de la siguiente manera y juzgar visualmente la gaussiness de los datos usando estas funciones (según la respuesta de VLC a la pregunta 2984)
Show[DensityHistogram[d, .2, ColorFunction -> (Opacity[0] &),
Method -> "DistributionAxes" -> "Histogram"], ListPlot[d]]
o estos:
GraphicsColumn[
Table[SmoothDensityHistogram[d, ColorFunction -> "DarkRainbow",
Method -> "DistributionAxes" -> p, ImageSize -> 500,
BaseStyle -> FontFamily -> "Helvetica",
LabelStyle -> Bold], p, True, "SmoothHistogram", "BoxWhisker"]]