2014-09-23 50 views
6

我有兩個變量,我已使用matplotlib散點圖功能繪製。 enter image description here二維圖的1sigma的置信區域

I would like to show the 68% confidence region by highlighting it in the plot.我知道,以顯示它在一個直方圖,但我不知道該怎麼做了這樣的2D圖(X對y)。在我的情況下,x is Massy is Ngal Mstar+2

的什麼,我找這個樣子的一個例子形象:

在這裏,他們都表現出使用深藍和95%置信區間採用淺藍色的68%置信區間。

是否可以使用scipy.stats模塊之一來實現?

enter image description here

+1

如果您可以訪問它,以適應數據,[seaborn](http://web.stanford.edu/ 〜mwaskom/software/seaborn/index.html)正是你想要的所有內容和內建 - 查找'regplot'。 – Ajean 2014-09-23 19:50:07

+0

@Ajean我相信你已經向我展示了我需要的東西!讓我試試吧 – ThePredator 2014-09-24 08:22:36

+0

我已經在python中使用了R包來做這件事,但我也對一個更簡單的解決方案感興趣。 – Doug 2014-10-21 18:59:58

回答

0

所有首先感謝您@snake_charmer的回答,但我發現從scipy.optimize

解決在使用curve_fit問題的一個簡單的方法我適合使用curve_fit這給了我最適合的參數我的數據樣本。它也給我的是參數的估計協方差。對角線的對角線提供了參數估計的方差。爲了計算參數上的一個標準差,我們可以使用np.sqrt(np.diag(pcov))其中pcov是協方差矩陣。

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

以上是我用於數據的適合函數。

現在使用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

要繪製的區域兩條曲線之間,你可以使用pyplot.fill_between()

至於你的信心區域,我不知道你想實現什麼,所以我舉例說明具有同步置信帶,通過修改代碼:

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

輸出如下這樣的: enter image description here