Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Question about stats.confband vs confbandl #3

Open
jrcooper91 opened this issue Jul 6, 2018 · 0 comments
Open

Question about stats.confband vs confbandl #3

jrcooper91 opened this issue Jul 6, 2018 · 0 comments

Comments

@jrcooper91
Copy link

jrcooper91 commented Jul 6, 2018

I have been using BCES in conjunction with nmmn for weighted linear regression and bootstrapping. It has been working fantastically so far. I have a question on usage between confband and confbandl. The documentation mentioned confband is for linear, while confbandl is nonlinear, the latter of which is used in the BCES example and modified for a weighted LR based on the function. I have been testing both stats modules, and I was surprised to see such a difference even when using BCES as my LR in both for the same dataset. Am I using one of these incorrectly? I originally thought to use stats.confbandl, but it seems from the documentation that stats.confband is more appropriate.

Edit: I also did this test with a random Gaussian distribution dataset, where the confbandl gave CI that were 2 orders of mag larger than expected, and the confband ones were within %5 of expectations. That's why I ultimately decided on utilizing confband. https://github.com/jrcooper91/van_der_Wel_size-mass/blob/master/Testing_bootstraps_BCES%20.ipynb
stats.confband

x = np.array([  0.0630552 ,0.138341,0.150978,1.509820,2.83197654,6.79663764,10.445637])
y = np.array([ 1.48,  1.14,  1.16,  1.15,  1.43,  1.43,  1.22])
yer = np.array([ 1.395,  0.1  ,  0.275,  0.12 ,  0.235,  0.135,  0.155])
sort = np.argsort(x)
x = x[sort]
x = np.array(x)
y = y[sort] 
y = np.array(y)
yer = yer[sort]
yer = np.array(yer)
xer=zeros(len(x))
cov=zeros(len(x))   # no correlation between error measurements
i=0 
nboot = 10000
a,b,erra,errb,covab=bces.bces.bcesp(x,xer,y,yer,cov,nboot)
ybces=a[3]*x+b[3]  
    # Gets lower and upper bounds on the confidence band 
lcb1,ucb1,x=nmmn.stats.confband(x, y, a[i], b[i], 0.68, x)
lcb2,ucb2,x2=nmmn.stats.confband(x, y, a[i], b[i], 0.95, x)
lcb3,ucb3,x3=nmmn.stats.confband(x, y, a[i], b[i], 0.997, x)
errorbar(x,y,yerr=yer,fmt='o')
ax = plot(x,ybces,'-k')
ax = fill_between(x, lcb1, ucb1, alpha=0.6, facecolor='purple')
ax = fill_between(x, lcb2, ucb2, alpha=0.3, facecolor='blue')
ax = fill_between(x, lcb3, ucb3, alpha=0.4, facecolor='grey')
ax = xlabel('x')
ax = ylabel('y')
plt.show()

statsconfband
stats.confbandl

x = np.array([ 0.0630552 , 0.13834198, 0.15097889, 1.50982026, 2.83197654, 6.79663764, 10.44563709])
y = np.array([ 1.48, 1.14, 1.16, 1.15, 1.43, 1.43, 1.22])
yer = np.array([ 1.395, 0.1 , 0.275, 0.12 , 0.235, 0.135, 0.155])
sort = np.argsort(x)
x = x[sort]
x = np.array(x)
y = y[sort]
y = np.array(y)
yer = yer[sort]
yer = np.array(yer)
xer=zeros(len(x))
cov=zeros(len(x)) # no correlation between error measurements
i=0
nboot = 10000
def func(x): return x[1]*x[0]+x[2]
a,b,erra,errb,covab=bces.bces.bcesp(x,xer,y,yer,cov,nboot)
ybces=a[3]*x+b[3] # the integer corresponds to the desired BCES method for plotting (3-ort, 0-y|x, 1-x|y, don't use bissector)
# array with best-fit parameters
fitm=np.array([ a[i],b[i] ])
# covariance matrix of parameter uncertainties
covm=np.array([ (erra[i]**2,covab[i]), (covab[i],errb[i]**2) ])
# Gets lower and upper bounds on the confidence band
lcb1,ucb1,xcb1=nmmn.stats.confbandnl(x,y,func,fitm,covm,2,0.68,x)
lcb2,ucb2,xcb2=nmmn.stats.confbandnl(x,y,func,fitm,covm,2,0.95, x)
lcb3,ucb3,xcb3=nmmn.stats.confbandnl(x,y,func,fitm,covm,2,0.997, x)
errorbar(xcb1,y,yerr=yer,fmt='o')
ax = plot(x,ybces,'-k')
ax = fill_between(x, lcb1, ucb1, alpha=0.6, facecolor='purple')
ax = fill_between(x, lcb2, ucb2, alpha=0.3, facecolor='blue')
ax = fill_between(x, lcb3, ucb3, alpha=0.4, facecolor='grey')
ax = xlabel('x')
ax = ylabel('y')
plt.show()

confbandl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant