source: orange/orange/Orange/regression/linear.py @ 7760:8f52f1b64337

Revision 7760:8f52f1b64337, 9.5 KB checked in by lanz <lan.zagar@…>, 3 years ago (diff)

Replaced function printLinearRegression with method str.

Line 
1from string import join
2
3import numpy
4from numpy.linalg import inv
5from numpy import dot, sqrt
6
7import Orange
8import statc
9
10
11########################################################################
12# Linear Regression
13
14class LinearRegressionLearner(object):
15    def __new__(self, data=None, name='linear regression', **kwds):
16        learner = object.__new__(self, **kwds)
17        if data:
18            learner.__init__(name,**kwds) # force init
19            return learner(data)
20        else:
21            return learner  # invokes the __init__
22
23    def __init__(self, name='linear regression', beta0=True,
24                 use_attributes=None, stepwise=False, add_sig=0.05,
25                 remove_sig=0.2, stepwise_before=True, **kw):
26        self.name = name
27        self.beta0 = beta0
28        self.stepwise = stepwise
29        self.stepwise_before = stepwise_before
30        self.add_sig = add_sig
31        self.remove_sig = remove_sig
32        self.use_attributes = use_attributes
33        self.__dict__.update(kw)
34
35    def __call__(self, data, weight=None):
36        if not self.use_attributes == None:
37            new_domain = Orange.data.Domain(self.use_attributes,
38                                            data.domain.classVar)
39            new_domain.addmetas(data.domain.getmetas())
40            data = Orange.data.Table(new_domain, data)
41           
42        if self.stepwise and self.stepwise_before:
43            use_attributes=stepwise(data, add_sig=self.add_sig,
44                                    remove_sig=self.remove_sig)
45            new_domain = Orange.data.Domain(use_attributes, data.domain.classVar)
46            new_domain.addmetas(data.domain.getmetas())
47            data = Orange.data.Table(new_domain, data)
48
49        # continuization (replaces discrete with continuous attributes)
50        continuizer = Orange.feature.continuization.DomainContinuizer()
51        continuizer.multinomialTreatment = continuizer.FrequentIsBase
52        continuizer.zeroBased = True
53        domain0 = continuizer(data)
54        data = data.translate(domain0)
55
56        if self.stepwise and not self.stepwise_before:
57            use_attributes = stepwise(data, weight, add_sig=self.add_sig,
58                                      remove_sig=self.remove_sig)
59            new_domain = Orange.data.Domain(use_attributes, data.domain.classVar)
60            new_domain.addmetas(data.domain.getmetas())
61            data = Orange.data.Table(new_domain, data)       
62       
63        # missing values handling (impute missing)
64        imputer = Orange.feature.imputation.ImputerConstructor_model()
65        imputer.learnerContinuous = Orange.regression.mean.MeanLearner()
66        imputer.learnerDiscrete = Orange.classification.majority.MajorityLearner()
67        imputer = imputer(data)
68        data = imputer(data)
69
70        # convertion to numpy
71        A, y, w = data.toNumpy()        # weights ??
72        if A == None:
73            n = len(data)
74            m = 0
75        else:
76            n, m = numpy.shape(A)
77     
78        if self.beta0 == True:
79             if A == None:
80                 X = numpy.ones([len(data), 1])
81             else:
82                 X = numpy.insert(A, 0, 1, axis=1) # adds a column of ones
83        else:
84             X = A
85
86        # set weights
87        W = numpy.identity(len(data))
88        if weight:
89            for di, d in enumerate(data):
90                W[di, di] = float(d[weight])
91
92        # adds some robustness by computing the pseudo inverse;
93        # normal inverse could fail due to singularity of the X.T * W * X
94        D = dot(dot(numpy.linalg.pinv(dot(dot(X.T, W), X)), X.T), W)
95        beta = dot(D, y)
96
97        yEstimated = dot(X, beta)  # estimation
98        # some desriptive statistisc
99        muY, sigmaY = numpy.mean(y), numpy.std(y)
100        muX, covX = numpy.mean(X, axis=0), numpy.cov(X, rowvar=0)
101 
102        # model statistics
103        SST = numpy.sum((y - muY) ** 2)
104        SSR = numpy.sum((yEstimated - muY) ** 2)
105        SSE, RSquare = SST - SSR, SSR / SST
106        R = numpy.sqrt(RSquare) # coefficient of determination
107        RAdjusted = 1 - (1 - RSquare) * (n - 1) / (n - m - 1)
108        F = (SSR / m) / (SST - SSR / (n - m - 1)) # F statistisc
109        df = m - 1
110 
111        sigmaSquare = SSE / (n - m - 1)
112
113        # standard error of estimated coefficients
114        errCoeff = sqrt(sigmaSquare * inv(dot(X.T, X)).diagonal())
115 
116        # t statistisc, significance
117        t = beta / errCoeff
118        df = n - 2
119        significance = []
120        for tt in t:
121            try:
122                significance.append(statc.betai(df * 0.5, 0.5,
123                                                df / (df + tt * tt)))
124            except:
125                significance.append(1.0)
126 
127        # standardized coefficients
128        if m > 0:
129            stdCoeff = (sqrt(covX.diagonal()) / sigmaY)  * beta
130        else:
131            stdCoeff = (sqrt(covX) / sigmaY)  * beta
132 
133        model = {'descriptives': { 'meanX': muX, 'covX': covX, 'meanY': muY,
134                                   'sigmaY': sigmaY},
135                 'model' : {'estCoeff': beta, 'stdErrorEstimation': errCoeff},
136                 'model summary': {'TotalVar': SST, 'ExplVar': SSE,
137                                   'ResVar': SSR, 'R': R, 'RAdjusted': RAdjusted,
138                                   'F': F, 't': t, 'sig': significance}}
139        return LinearRegression(statistics=model, domain=data.domain,
140                                name=self.name, beta0=self.beta0, imputer=imputer)
141
142class LinearRegression(Orange.classification.Classifier):
143    def __init__(self, **kwds):
144        for a, b in kwds.items():
145            self.setattr(a, b)
146        self.beta = self.statistics['model']['estCoeff']
147
148    def __call__(self, example, resultType=Orange.classification.Classifier.GetValue):
149        ex = Orange.data.Instance(self.domain, example)
150        ex = self.imputer(ex)
151        ex = numpy.array(ex.native())
152
153        if self.beta0:
154            if len(self.beta) > 1:
155                yhat = self.beta[0] + dot(self.beta[1:], ex[:-1])
156            else:
157                yhat = self.beta[0]
158        else:
159            yhat = dot(self.beta, ex[:-1])
160        yhat = Orange.data.Value(yhat)
161         
162        if resultType == Orange.classification.Classifier.GetValue:
163            return yhat
164        if resultType == Orange.classification.Classifier.GetProbabilities:
165            return Orange.statistics.distribution.Continuous({1.0: yhat})
166        return (yhat, Orange.statistics.distribution.Continuous({1.0: yhat}))
167   
168    def __str__(self):
169        err = self.statistics['model']['stdErrorEstimation']
170        t = self.statistics['model summary']['t']
171        sig = self.statistics['model summary']['sig']
172       
173        s = ' '.join(['%10s' % l for l in
174                      ('Variable', 'Coeff Est', 'Std Error', 't-value', 'p')])
175
176        fmt = '\n%10s ' + ' '.join(["%10.3f"] * 4)
177        if self.beta0 == True:
178            s += fmt % ('Constant', self.beta[0], err[0], t[0], sig[0])
179            for i in range(len(self.domain.attributes) - 1):
180                s +=  fmt % (self.domain.attributes[i].name,
181                             self.beta[i + 1], err[i + 1], t[i + 1], sig[i + 1])
182        else:
183            for i in range(len(self.domain.attributes) - 1):
184                s +=  fmt % (self.domain.attributes[i].name,
185                             self.beta[i], err[i], t[i], sig[i])
186        return s
187
188def get_sig(m1, m2, n):
189    if m1 == None or m2 == None:
190        return 1.0
191    p1, p2 = len(m1.domain.attributes), len(m2.domain.attributes)
192    RSS1 = m1.statistics["model summary"]["ExplVar"]
193    RSS2 = m2.statistics["model summary"]["ExplVar"]
194    if RSS1 <= RSS2 or p2 <= p1 or n <= p2 or RSS2 <= 0:
195        return 1.0
196    F = ((RSS1 - RSS2) / (p2 - p1)) / (RSS2 / (n - p2))
197    return statc.fprob(int(p2 - p1), int(n - p2), F)
198
199def stepwise(data, weight, add_sig=0.05, remove_sig=0.2):
200    inc_atts = []
201    not_inc_atts = [at for at in data.domain.attributes]
202
203    changed_model = True
204    while changed_model:
205        changed_model = False
206        # remove all unsignificant conditions (learn several models,
207        # where each time one attribute is removed and check significance)
208        orig_lin_reg = LinearRegressionLearner(data, use_attributes=inc_atts)
209        reduced_lin_reg = []
210        for ati in range(len(inc_atts)):
211            try:
212                reduced_lin_reg.append(LinearRegressionLearner(data, weight,
213                        use_attributes=inc_atts[:ati] + inc_atts[(ati + 1):]))
214            except:
215                reduced_lin_reg.append(None)
216       
217        sigs = [get_sig(r, orig_lin_reg, len(data)) for r in reduced_lin_reg]
218        if sigs and max(sigs) > remove_sig:
219            # remove that attribute, start again
220            crit_att = inc_atts[sigs.index(max(sigs))]
221            not_inc_atts.append(crit_att)
222            inc_atts.remove(crit_att)
223            changed_model = True
224            continue
225
226        # add all significant conditions (go through all attributes in
227        # not_inc_atts, is there one that significantly improves the model?
228        more_complex_lin_reg = []
229        for ati in range(len(not_inc_atts)):
230            try:
231                more_complex_lin_reg.append(LinearRegressionLearner(data,
232                        weight, use_attributes=inc_atts + [not_inc_atts[ati]]))
233            except:
234                more_complex_lin_reg.append(None)
235
236        sigs = [get_sig(orig_lin_reg, r, len(data)) for r in more_complex_lin_reg]
237        if sigs and min(sigs) < add_sig:
238            best_att = not_inc_atts[sigs.index(min(sigs))]
239            inc_atts.append(best_att)
240            not_inc_atts.remove(best_att)
241            changed_model = True
242    return inc_atts
243
Note: See TracBrowser for help on using the repository browser.