1sigma trust areas for 2D graphics

I have two variables that I built using the matplotlib scatter matplotlib . enter image description here

I would like to show the 68% confidence region by highlighting it in the plot. I know to show it in a histogram, but I don't know how to do it for a 2D plot like this (x vs y). In my case, x is Mass and y is Ngal Mstar+2 .

An example image of what I'm looking for is as follows:

Here they showed a 68% confidence area using navy blue and a 95% confidence range using light blue.

Can this be achieved using one of the scipy.stats modules ?

enter image description here

+6
source share
2 answers

First of all, thanks @snake_charmer for your answer , but I found an easier way to solve the problem using curve_fit from scipy.optimize

I customize my sample data using curve_fit , which gives me my best suitable parameters. What this also gives me is an estimate of the covariance of the parameters. Diagonals of the same parameter provide variance of the parameter estimate. To calculate one standard error of the deviation in the parameters, we can use np.sqrt(np.diag(pcov)) , where pcov is the covariance matrix.

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

The above fit function that I use for data.

Now to put data using 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

0
source

To plot between two curves, you can use pyplot.fill_between() .

As for your region of trust, I was not sure what you wanted to achieve, so I demonstrated simultaneous confidence ranges by changing the code:

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()/(nX.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], nX.shape[1], 0.95) QF_2 = sp.fdtri(X.shape[1], nX.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() 

The result is as follows: enter image description here

+4
source

Source: https://habr.com/ru/post/975709/


All Articles