Orchid XP v8 wrote:
>>> Any idea how to get a confidence estimate?
>>
>> You could take a look at the section labeled "Confidence Intervals" on 
>> the page that scott linked.
> 
> Yes. Except, as originally explained, I don't understand it.
> 

Considering that you're only looking for an implementation it doesn't 
really matter that you understand what's going on so long as you can 
just type in the equations.  Everything is pretty straightforward except 
for determining t*_(n-2), for which I suggest either finding an existing 
implementation of the inverse cumulative Student's t-distribution 
function or looking up the appropriate value from the table on 
http://en.wikipedia.org/wiki/Student's_t-distribution.

In case you're still having trouble, I whipped up some Python code to do 
what you want.  It's almost pseudocode, so hopefully it'll be pretty 
easy to follow even if you don't know Python.


------

from math import sqrt
from scipy.stats import t as tFcn


def linRegress(X, Y, conf):
	# maximum likelihood estimates
	n = len(X)
	Sx = sum(X)
	Sy = sum(Y)
	Sxx = sum(x*x for x in X)
	Sxy = sum(x*y for x, y in zip(X,Y))
	Syy = sum(y*y for y in Y)
	b = (n*Sxy - Sx*Sy)/(n*Sxx - Sx*Sx)
	a = Sy/n - b*Sx/n

	# confidence interval stuff
	ve = (n*Syy - Sy*Sy - b*b*(n*Sxx - Sx*Sx))/(n*(n-2))
	vb = n*ve / (n*Sxx - Sx*Sx)
	va = vb*Sxx / n

	tstar = tFcn(n-2).isf(1-conf)
	return (a, tstar*sqrt(va)), (b, tstar*sqrt(vb))


if __name__ == "__main__":
	# example data from the Wikipedia page
	X = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 
1.75, 1.78, 1.80, 1.83]
	Y = [52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 
64.47, 66.28, 68.10, 69.92, 72.19, 74.46]

	a, b = linRegress(X,Y,0.975)
	print "intercept:  mean=%f min=%f max=%f" % (a[0], a[0]-a[1], a[0]+a[1])
	print "slope:  mean=%f min=%f max=%f" % (b[0], b[0]-b[1], b[0]+b[1])

Post a reply to this message