Calculating the prediction band of a linear regression model

The method below calculates the prediction band of an arbitrary linear regression model at a given confidence level in Python.

If you use it, let me know if you find any bugs.

Note: the scatterfit method called below is available here.

2. def predband(xd,yd,a,b,conf=0.95,x=None):
3.     """
4. Calculates the prediction band of the linear regression model at the desired confidence
5. level.
7. Clarification of the difference between confidence and prediction bands:
8. "The 2sigma confidence interval is 95% sure to contain the best-fit regression line.
9. This is not the same as saying it will contain 95% of the data points. The prediction bands are
10. further from the best-fit line than the confidence bands, a lot further if you have many data
11. points. The 95% prediction interval is the area in which you expect 95% of all data points to fall."
12. (from http://graphpad.com/curvefit/linear_regression.htm)
14. Arguments:
15. - conf: desired confidence level, by default 0.95 (2 sigma)
16. - xd,yd: data arrays
17. - a,b: linear fit parameters as in y=ax+b
18. - x: (optional) array with x values to calculate the confidence band. If none is provided, will
19.  by default generate 100 points in the original x-range of the data.
20.  
21. Usage:
22. >>> lpb,upb,x=nemmen.predband(all.kp,all.lg,a,b,conf=0.95)
23. calculates the prediction bands for the given input arrays
25. >>> pylab.fill_between(x, lpb, upb, alpha=0.3, facecolor='gray')
26. plots a shaded area containing the prediction band  
28. Returns:
29. Sequence (lpb,upb,x) with the arrays holding the lower and upper confidence bands
30. corresponding to the [input] x array.
32. References:
33. 1. http://www.JerryDallal.com/LHSP/slr.htm, Introduction to Simple Linear Regression, Gerard
34. E. Dallal, Ph.D.
36. Rodrigo Nemmen
37. v1 Dec. 2011
38. v2 Jun. 2012: corrected bug in dy.
39.     """
40.     alpha=1.-conf   # significance
41.     n=xd.size   # data sample size
43.     if x==None: x=numpy.linspace(xd.min(),xd.max(),100)
45.     # Predicted values (best-fit model)
46.     y=a*x+b
48.     # Auxiliary definitions
49.     sd=scatterfit(xd,yd,a,b)    # Scatter of data about the model
50.     sxd=numpy.sum((xd-xd.mean())**2)
51.     sx=(x-xd.mean())**2 # array
53.     # Quantile of Student's t distribution for p=1-alpha/2
54.     q=scipy.stats.t.ppf(1.-alpha/2.,n-2)
56.     # Prediction band
57.     dy=q*sd*numpy.sqrt( 1.+1./n + sx/sxd )
58.     upb=y+dy    # Upper prediction band
59.     lpb=y-dy    # Lower prediction band
61.     return lpb,upb,x