2014-09-23 11 views
6

Ho due variabili che ho tracciato usando la funzione scatter matplotlib. enter image description hereAree di confidenza di 1sigma per un grafico 2D

I would like to show the 68% confidence region by highlighting it in the plot. So di mostrarlo in un istogramma, ma non so come si fa per una trama 2D come questo (x vs y). Nel mio caso, il x is Mass e y is Ngal Mstar+2.

Un'immagine esempio di quello che sto cercando per un look come questo:

Qui hanno mostrato la regione di confidenza 68% utilizzando blu e il 95% scuro regione di confidenza con la luce blu.

Può essere ottenuto utilizzando uno dei moduliscipy.stats?

enter image description here

+1

Se si ha accesso ad esso, [Seaborn] (http://web.stanford.edu/ ~ mwaskom/software/seaborn/index.html) ha proprio tutto ciò che desideri e tutto ciò che vuoi - guarda su 'regplot'. – Ajean

+0

@Ajean Credo che tu mi abbia mostrato la cosa di cui ho bisogno! Fammi provare – ThePredator

+0

Ho usato pacchetti R in python per fare questo in passato, ma mi interessa anche una soluzione più semplice. – Doug

risposta

0

Prima di tutto grazie @snake_charmer per la risposta, ma ho trovato un modo più semplice per risolvere il problema utilizzandocurve_fit da scipy.optimize

mi adatto il mio campione di dati utilizzando curve_fit che mi dà i miei migliori parametri di fit. Ciò che mi dà anche è la stima covarianza dei parametri. Le diagonali della stessa forniscono la varianza della stima dei parametri. Per calcolare uno degli errori di deviazione standard sui parametri, è possibile utilizzare np.sqrt(np.diag(pcov)) dove pcov è la matrice di covarianza.

def fitfunc(M,p1,p2): 
    N = p1+((M)*p2) 
    return N 

Quanto sopra è la funzione di adattamento che utilizzo per i dati.

Ora per adattarsi ai dati utilizzando curve_fit

popt_1,pcov_1 = curve_fit(fitfunc,logx,logn,p0=(10.0,1.0),maxfev=2000) 

p1_1 = popt_1[0] 
p1_2 = popt_1[1] 

sigma1 = [np.sqrt(pcov_1[0,0]),np.sqrt(pcov_1[1,1])] #THE 1 SIGMA CONFIDENCE INTERVALS 
residuals1 = (logy) - fitfunc((logx),p1_1,p1_2) 
xi_sq_1 = sum(residuals1**2) #THE CHI-SQUARE OF THE FIT 

curve_y_1 = fitfunc((logx),p1_1,p1_2) 

fig = plt.figure() 
ax1 = fig.add_subplot(111) 
ax1.scatter(logx,logy,c='r',label='$0.0<z<0.5$') 
ax1.plot(logx,curve_y_1,'y') 
ax1.plot(logx,fitfunc(logx,p1_1+sigma1[0],p1_2+sigma1[1]),'m',label='68% conf limits') 
ax1.plot(logx,fitfunc(logx,p1_1-sigma1[0],p1_2-sigma1[1]),'m') 

So just by using the square root the diagonal elements of the covariance matrix, I can obtain the 1 sigma confidence lines.

enter image description here

4

per tracciare una regione tra due curve, è possibile utilizzare pyplot.fill_between().

quanto riguarda la tua regione di confidenza, non ero sicuro di quello che si voleva ottenere, così ho esemplificato con bande di confidenza simultanei, modificando il codice da:

https://en.wikipedia.org/wiki/Confidence_and_prediction_bands#cite_note-2

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.special as sp 

## Sample size. 
n = 50 

## Predictor values. 
XV = np.random.uniform(low=-4, high=4, size=n) 
XV.sort() 

## Design matrix. 
X = np.ones((n,2)) 
X[:,1] = XV 

## True coefficients. 
beta = np.array([0, 1.], dtype=np.float64) 

## True response values. 
EY = np.dot(X, beta) 

## Observed response values. 
Y = EY + np.random.normal(size=n)*np.sqrt(20) 

## Get the coefficient estimates. 
u,s,vt = np.linalg.svd(X,0) 
v = np.transpose(vt) 
bhat = np.dot(v, np.dot(np.transpose(u), Y)/s) 

## The fitted values. 
Yhat = np.dot(X, bhat) 

## The MSE and RMSE. 
MSE = ((Y-EY)**2).sum()/(n-X.shape[1]) 
s = np.sqrt(MSE) 

## These multipliers are used in constructing the intervals. 
XtX = np.dot(np.transpose(X), X) 
V = [np.dot(X[i,:], np.linalg.solve(XtX, X[i,:])) for i in range(n)] 
V = np.array(V) 

## The F quantile used in constructing the Scheffe interval. 
QF = sp.fdtri(X.shape[1], n-X.shape[1], 0.95) 
QF_2 = sp.fdtri(X.shape[1], n-X.shape[1], 0.68) 

## The lower and upper bounds of the Scheffe band. 
D = s*np.sqrt(X.shape[1]*QF*V) 
LB,UB = Yhat-D,Yhat+D 
D_2 = s*np.sqrt(X.shape[1]*QF_2*V) 
LB_2,UB_2 = Yhat-D_2,Yhat+D_2 


## Make the plot. 
plt.clf() 
plt.plot(XV, Y, 'o', ms=3, color='grey') 
plt.hold(True) 
a = plt.plot(XV, EY, '-', color='black', zorder = 4) 

plt.fill_between(XV, LB_2, UB_2, where = UB_2 >= LB_2, facecolor='blue', alpha= 0.3, zorder = 0) 
b = plt.plot(XV, LB_2, '-', color='blue', zorder=1) 
plt.plot(XV, UB_2, '-', color='blue', zorder=1) 

plt.fill_between(XV, LB, UB, where = UB >= LB, facecolor='blue', alpha= 0.3, zorder = 2) 
b = plt.plot(XV, LB, '-', color='blue', zorder=3) 
plt.plot(XV, UB, '-', color='blue', zorder=3) 

d = plt.plot(XV, Yhat, '-', color='red',zorder=4) 

plt.ylim([-8,8]) 
plt.xlim([-4,4]) 

plt.xlabel("X") 
plt.ylabel("Y") 

plt.show() 

L'uscita sembra questo: enter image description here