intervalos de confianza multidimensionales

Tengo numerosas tuplas (par1, par2), es decir, puntos en un espacio de parámetros bidimensionales obtenidos al repetir un experimento varias veces.

Estoy buscando la posibilidad de calcular y visualizar las elipses de confianza (no estoy seguro si ese es el término correcto para esto). Aquí una ttwig de ejemplo que encontré en la web para mostrar lo que quiero decir:

introduzca la descripción de la imagen aquí

fuente: blogspot.ch/2011/07/classification-and-discrimination-with.html

Entonces, en principio, uno tiene que ajustar una distribución normal multivariable a un histogtwig 2D de puntos de datos, supongo. puede alguien ayudarme con esto?

Parece que solo quieres la elipse 2-sigma de la dispersión de puntos?

Si es así, considere algo como esto (del código para un documento aquí: https://github.com/joferkington/oost_paper_code/blob/master/error_ellipse.py ):

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse def plot_point_cov(points, nstd=2, ax=None, **kwargs): """ Plots an `nstd` sigma ellipse based on the mean and covariance of a point "cloud" (points, an Nx2 array). Parameters ---------- points : An Nx2 array of the data points. nstd : The radius of the ellipse in numbers of standard deviations. Defaults to 2 standard deviations. ax : The axis that the ellipse will be plotted on. Defaults to the current axis. Additional keyword arguments are pass on to the ellipse patch. Returns ------- A matplotlib ellipse artist """ pos = points.mean(axis=0) cov = np.cov(points, rowvar=False) return plot_cov_ellipse(cov, pos, nstd, ax, **kwargs) def plot_cov_ellipse(cov, pos, nstd=2, ax=None, **kwargs): """ Plots an `nstd` sigma error ellipse based on the specified covariance matrix (`cov`). Additional keyword arguments are passed on to the ellipse patch artist. Parameters ---------- cov : The 2x2 covariance matrix to base the ellipse on pos : The location of the center of the ellipse. Expects a 2-element sequence of [x0, y0]. nstd : The radius of the ellipse in numbers of standard deviations. Defaults to 2 standard deviations. ax : The axis that the ellipse will be plotted on. Defaults to the current axis. Additional keyword arguments are pass on to the ellipse patch. Returns ------- A matplotlib ellipse artist """ def eigsorted(cov): vals, vecs = np.linalg.eigh(cov) order = vals.argsort()[::-1] return vals[order], vecs[:,order] if ax is None: ax = plt.gca() vals, vecs = eigsorted(cov) theta = np.degrees(np.arctan2(*vecs[:,0][::-1])) # Width and height are "full" widths, not radius width, height = 2 * nstd * np.sqrt(vals) ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, **kwargs) ax.add_artist(ellip) return ellip if __name__ == '__main__': #-- Example usage ----------------------- # Generate some random, correlated data points = np.random.multivariate_normal( mean=(1,1), cov=[[0.4, 9],[9, 10]], size=1000 ) # Plot the raw points... x, y = points.T plt.plot(x, y, 'ro') # Plot a transparent 3 standard deviation covariance ellipse plot_point_cov(points, nstd=3, alpha=0.5, color='green') plt.show() 

introduzca la descripción de la imagen aquí

Consulte la publicación Cómo dibujar una elipse de error de covarianza .

Aquí está la realización de python:

 import numpy as np from scipy.stats import norm, chi2 def cov_ellipse(cov, q=None, nsig=None, **kwargs): """ Parameters ---------- cov : (2, 2) array Covariance matrix. q : float, optional Confidence level, should be in (0, 1) nsig : int, optional Confidence level in unit of standard deviations. Eg 1 stands for 68.3% and 2 stands for 95.4%. Returns ------- width, height, rotation : The lengths of two axises and the rotation angle in degree for the ellipse. """ if q is not None: q = np.asarray(q) elif nsig is not None: q = 2 * norm.cdf(nsig) - 1 else: raise ValueError('One of `q` and `nsig` should be specified.') r2 = chi2.ppf(q, 2) val, vec = np.linalg.eigh(cov) width, height = 2 * sqrt(val[:, None] * r2) rotation = np.degrees(arctan2(*vec[::-1, 0])) return width, height, rotation 

El significado de la desviación estándar es incorrecto en la respuesta de Joe Kington. Usualmente usamos 1, 2 sigma para 68%, 95% de niveles de confianza, pero la elipse de 2 sigma en su respuesta no contiene 95% de probabilidad de la distribución total. La forma correcta es usar una distribución de chi cuadrado para estimar el tamaño de la elipse como se muestra en la publicación .

Modifiqué ligeramente uno de los ejemplos anteriores que traza los contornos de la región de confianza o de error. Ahora creo que le da los contornos correctos.

Estaba dando los contornos incorrectos porque estaba aplicando el método scoreatpercentile al conjunto de datos conjunto (puntos azules + rojos) cuando debería aplicarse por separado a cada conjunto de datos.

El código modificado se puede encontrar a continuación:

 import numpy import scipy import scipy.stats import matplotlib.pyplot as plt # generate two normally distributed 2d arrays x1=numpy.random.multivariate_normal((100,420),[[120,80],[80,80]],400) x2=numpy.random.multivariate_normal((140,340),[[90,-70],[-70,80]],400) # fit a KDE to the data pdf1=scipy.stats.kde.gaussian_kde(x1.T) pdf2=scipy.stats.kde.gaussian_kde(x2.T) # create a grid over which we can evaluate pdf q,w=numpy.meshgrid(range(50,200,10), range(300,500,10)) r1=pdf1([q.flatten(),w.flatten()]) r2=pdf2([q.flatten(),w.flatten()]) # sample the pdf and find the value at the 95th percentile s1=scipy.stats.scoreatpercentile(pdf1(pdf1.resample(1000)), 5) s2=scipy.stats.scoreatpercentile(pdf2(pdf2.resample(1000)), 5) # reshape back to 2d r1.shape=(20,15) r2.shape=(20,15) # plot the contour at the 95th percentile plt.contour(range(50,200,10), range(300,500,10), r1, [s1],colors='b') plt.contour(range(50,200,10), range(300,500,10), r2, [s2],colors='r') # scatter plot the two normal distributions plt.scatter(x1[:,0],x1[:,1],alpha=0.3) plt.scatter(x2[:,0],x2[:,1],c='r',alpha=0.3) 

Supongo que lo que está buscando es calcular las regiones de confianza .

No sé mucho sobre eso, pero como punto de partida, revisaría la aplicación del sherpa para Python. Al menos, en su charla sobre Scipy 2011, los autores mencionan que puede determinar y obtener regiones de confianza con ella (aunque es posible que necesite un modelo para sus datos).

Vea el video y las diapositivas correspondientes de la charla de Sherpa.

HTH