|
 |
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
|
|
|
|
Invisible wrote:
> We'll see if I can deduce anything further...
Public Function FACTLN(N)
Dim A(1 To 100) As Double
Dim I As Integer
For I = 1 To 100
A(I) = -100
Next I
If N < 0 Then End
If N <= 99 Then
If A(N + 1) < 0 Then A(N + 1) = GAMMALN(N + 1)
FACTLN = A(N + 1)
Else
FACTLN = GAMMALN(N + 1)
End If
End Function
Uh... anybody want to venture a guess is to how this is different from
Public Function FACTLN(N)
IF N < 0 Then End
FACTLN = GAMMALN(N + 1)
End Function
From what I can tell, it creates an array, initialises it with negative
values, then checks to see if a given cell is negative (which it will
be, 100% of the time), and if so, calculates a value, puts it into the
array, and returns it.
It _looks like_ it's trying to cache previously-computed results. Except
that the array gets recreated each time the function is called. (Unless
I've missed something.) And even if it doesn't, each call to the
function zeros out the whole array anyway! WTF?
Post a reply to this message
|
|
