source: orange/orange/Orange/regression/linear.py @ 7762:68723affedca

Revision 7762:68723affedca, 9.5 KB checked in by lanz <lan.zagar@…>, 3 years ago (diff)

Changed some names to use underscores.

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