|
![](/i/fill.gif) |
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
|
![](/i/fill.gif) |