1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
"""
General linear models
--------------------
"""
import numpy as N
from scipy.sandbox.models import family
from scipy.sandbox.models.regression import wls_model
class model(wls_model):
niter = 10
def __init__(self, design, family=family.Gaussian()):
self.family = family
super(model, self).__init__(design, weights=1)
def __iter__(self):
self.iter = 0
self.dev = N.inf
return self
def deviance(self, Y=None, results=None, scale = 1.):
"""
Return (unnormalized) log-likelihood for glm.
Note that self.scale is interpreted as a variance in old_model, so
we divide the residuals by its sqrt.
"""
if results is None:
results = self.results
if Y is None:
Y = self.Y
return self.family.deviance(Y, results.mu) / scale
def next(self):
results = self.results
Y = self.Y
self.weights = self.family.weights(results.mu)
self.initialize(self.design)
Z = results.predict + self.family.link.deriv(results.mu) * (Y - results.mu)
newresults = super(model, self).fit(self, Z)
newresults.Y = Y
newresults.mu = self.family.link.inverse(newresults.predict)
self.iter += 1
return newresults
def cont(self, tol=1.0e-05):
"""
Continue iterating, or has convergence been obtained?
"""
if self.iter >= model.niter:
return False
curdev = self.deviance(results=self.results)
if N.fabs((self.dev - curdev) / curdev) < tol:
return False
self.dev = curdev
return True
def estimate_scale(self, Y=None, results=None):
"""
Return Pearson\'s X^2 estimate of scale.
"""
if results is None:
results = self.results
if Y is None:
Y = self.Y
resid = Y - results.mu
return ((N.power(resid, 2) / self.family.variance(results.mu)).sum()
/ results.df_resid)
def fit(self, Y):
self.Y = N.asarray(Y, N.float64)
iter(self)
self.results = super(model, self).fit(
self.family.link.initialize(Y))
self.results.mu = self.family.link.inverse(self.results.predict)
self.scale = self.results.scale = self.estimate_scale()
while self.cont(self.results):
self.results = self.next()
self.scale = self.results.scale = self.estimate_scale()
return self.results
|
