Obszary ufności 1sigma dla wykresu 2D

Question

Mam dwie zmienne, które zostały wykreślone przy użyciu funkcji rozproszenia matplotlib.

I would like to show the 68% confidence region by highlighting it in the plot. Znam go pokazać w histogramie, ale nie wiem jak to zrobić za pomocą 2D działki tak (x vs y). W moim przypadku: x is Mass i y is Ngal Mstar+2.

Przykładowy obraz tego, co szukam wygląda następująco:

Oto one wykazały 68% obszaru ufności używając granatową i 95% ufności regionie za pomocą jasnoniebieskim.

Czy można to osiągnąć za pomocą jednego z modułówscipy.stats?

Answer 1

Przede wszystkim dziękuję @snake_charmer za odpowiedź, ale znalazłem prostszy sposób rozwiązania problemu za pomocącurve_fit z scipy.optimize

pasuję moje próbki danych przy użyciu curve_fit co daje mi moje najlepsze parametry pasują. To, co mi daje, to szacowana kowariancja parametrów. Przekątne tego samego zapewniają wariancję estymacji parametru. Aby obliczyć jeden błąd odchylenia standardowego na parametrach, można użyć np.sqrt(np.diag(pcov)), gdzie pcov jest macierzą kowariancji.

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

Powyższa funkcja służy do dopasowania danych.

Teraz, aby pasowały do ​​danych za pomocą 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.

Answer 2

Aby wykreślić obszar między dwiema krzywymi, można użyć pyplot.fill_between().

Jak dla swojego regionu ufności, nie byłem pewien, co chcesz osiągnąć, więc przykładowo z jednoczesnymi zespołów ufności, modyfikując kod:

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() 

Wyjście wygląda to: