# source:orange/Orange/orng/orngRegression.py@9671:a7b056375472

Revision 9671:a7b056375472, 12.6 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Moved orange to Orange (part 2)

Line
1import orange
2import numpy
3from numpy.linalg import inv
4from numpy import dot, sqrt
5import statc
6from string import join
7
8########################################################################
9# Linear Regression
10
11class LinearRegressionLearner(object):
12    def __new__(self, data=None, name='linear regression', **kwds):
13        learner = object.__new__(self, **kwds)
14        if data:
15            learner.__init__(name,**kwds) # force init
16            return learner(data)
17        else:
18            return learner  # invokes the __init__
19
20    def __init__(self, name='linear regression', beta0 = True, use_attributes = None, stepwise=False, add_sig=0.05, remove_sig=0.2, stepwise_before = True, **kw):
21        self.name = name
22        self.beta0 = beta0
23        self.stepwise=stepwise
24        self.stepwise_before = stepwise_before
26        self.remove_sig=remove_sig
27        self.use_attributes = use_attributes
28        self.__dict__.update(kw)
29
30    def __call__(self, data, weight=None):
31        if not self.use_attributes == None:
32            new_domain = orange.Domain(self.use_attributes, data.domain.classVar)
34            data = orange.ExampleTable(new_domain, data)
35
36        if self.stepwise and self.stepwise_before:
38            new_domain = orange.Domain(use_attributes, data.domain.classVar)
40            data = orange.ExampleTable(new_domain, data)
41
42        # continuization (replaces discrete with continuous attributes)
43        continuizer = orange.DomainContinuizer()
44        continuizer.multinomialTreatment = continuizer.FrequentIsBase
45        continuizer.zeroBased = True
46        domain0 = continuizer(data)
47        data = data.translate(domain0)
48
49        if self.stepwise and not self.stepwise_before:
51            new_domain = orange.Domain(use_attributes, data.domain.classVar)
53            data = orange.ExampleTable(new_domain, data)
54
55        # missing values handling (impute missing)
56        imputer = orange.ImputerConstructor_model()
57        imputer.learnerContinuous = orange.MajorityLearner()
58        imputer.learnerDiscrete = orange.MajorityLearner()
59        imputer = imputer(data)
60        data = imputer(data)
61
62        # convertion to numpy
63        A, y, w = data.toNumpy()        # weights ??
64        if A==None:
65            n = len(data)
66            m = 0
67        else:
68            n, m = numpy.shape(A)
69
70        if self.beta0 == True:
71             if A==None:
72                 X = numpy.ones([len(data),1])
73             else:
74                 X = numpy.insert(A,0,1,axis=1) # adds a column of ones
75        else:
76             X = A
77
78        # set weights
79        W = numpy.identity(len(data))
80        if weight:
81            for di, d in enumerate(data):
82                W[di,di] = float(d[weight])
83
84        D = dot(dot(numpy.linalg.pinv(dot(dot(X.T,W),X)), X.T), W) # adds some robustness by computing the pseudo inverse; normal inverse could fail due to singularity of the X.T*W*X
85        beta = dot(D,y)
86
87        yEstimated = dot(X,beta)  # estimation
88        # some desriptive statistisc
89        muY, sigmaY = numpy.mean(y), numpy.std(y)
90        muX, covX = numpy.mean(X, axis = 0), numpy.cov(X, rowvar = 0)
91
92        # model statistics
93        SST, SSR = numpy.sum((y - muY) ** 2), numpy.sum((yEstimated - muY) ** 2)
94        SSE, RSquare = SST-SSR, SSR/SST
95        R = numpy.sqrt(RSquare) # coefficient of determination
96        RAdjusted = 1 - (1 - RSquare) * (n - 1) / (n - m - 1)
97        F = (SSR / m) / (SST - SSR / (n - m - 1)) # F statistisc
98        df = m - 1
99
100        sigmaSquare = SSE / (n-m-1)
101
102        # standard error of estimated coefficients
103        errCoeff = sqrt(sigmaSquare * inv(dot(X.T,X)).diagonal())
104
105        # t statistisc, significance
106        t = beta / errCoeff
107        df = n-2
108        significance = []
109        for tt in t:
110            try:
111                significance.append(statc.betai(df*0.5,0.5,df/(df+tt*tt)))
112            except:
113                significance.append(1.0)
114
115        # standardized coefficients
116        if m>0:
117            stdCoeff = (sqrt(covX.diagonal()) / sigmaY)  * beta
118        else:
119            stdCoeff = (sqrt(covX) / sigmaY)  * beta
120
121        model = {'descriptives': { 'meanX' : muX, 'covX' : covX, 'meanY' : muY, 'sigmaY' : sigmaY},
122                 'model' : {'estCoeff' : beta, 'stdErrorEstimation': errCoeff},
123                 'model summary': {'TotalVar' : SST, 'ExplVar' : SSE, 'ResVar' : SSR, 'R' : R, 'RAdjusted' : RAdjusted,
124                                   'F' : F, 't' : t, 'sig': significance}}
125        return LinearRegression(statistics = model, domain = data.domain, name = self.name, beta0 = self.beta0, imputer=imputer)
126
127class LinearRegression(orange.Classifier):
128    def __init__(self, **kwds):
129        for a,b in kwds.items():
130            self.setattr(a,b)
131        self.beta = self.statistics['model']['estCoeff']
132
133    def __call__(self, example, resultType = orange.GetValue):
134        ex = orange.Example(self.domain, example)
135        ex = self.imputer(ex)
136        ex = numpy.array(ex.native())
137
138        if self.beta0:
139            if len(self.beta)>1:
140                yhat = self.beta[0] + dot(self.beta[1:], ex[:-1])
141            else:
142                yhat = self.beta[0]
143        else:
144            yhat = dot(self.beta, ex[:-1])
145        yhat = orange.Value(yhat)
146
147        if resultType == orange.GetValue:
148            return yhat
149        if resultType == orange.GetProbabilities:
150            return orange.ContDistribution({1.0: yhat})
151        return (yhat, orange.ContDistribution({1.0: yhat}))
152
153
154def printLinearRegression(lr):
155    """pretty-prints linear regression model"""
156    beta = lr.beta
157    err = lr.statistics['model']['stdErrorEstimation']
158    t = lr.statistics['model summary']['t']
159    sig = lr.statistics['model summary']['sig']
160    beta0 = lr.beta0
161
162    labels = ('Variable', 'Coeff Est', 'Std Error', 't-value', 'p')
163    print join(['%10s' % l for l in labels], ' ')
164
165    fmt = "%10s " + join(["%10.3f"]*4, " ")
166    if beta0 == True:
167        print fmt % ('Constant', beta[0], err[0], t[0], sig[0])
168        for i in range(len(lr.domain.attributes)-1):
169            print fmt % (lr.domain.attributes[i].name, beta[i+1], err[i+1], t[i+1], sig[i+1])
170    else:
171        for i in range(len(lr.domain.attributes)-1):
172            print fmt % (lr.domain.attributes[i].name, beta[i], err[i], t[i], sig[i])
173
174def get_sig(m1, m2, n):
175    if m1==None or m2==None:
176        return 1.0
177    p1, p2 = len(m1.domain.attributes), len(m2.domain.attributes)
180        return 1.0
182    return statc.fprob(int(p2-p1),int(n-p2),F)
183
184def stepwise(data, weight, add_sig = 0.05, remove_sig = 0.2):
185    inc_atts = []
186    not_inc_atts = [at for at in data.domain.attributes]
187
188    changed_model = True
189    while changed_model:
190        changed_model = False
191        # remove all unsignificant conditions (learn several models, where each time
192        # one attribute is removed and check significance)
193        orig_lin_reg = LinearRegressionLearner(data, use_attributes = inc_atts)
194        reduced_lin_reg = []
195        for ati in range(len(inc_atts)):
196            try:
197                reduced_lin_reg.append(LinearRegressionLearner(data, weight, use_attributes = inc_atts[:ati]+inc_atts[(ati+1):]))
198            except:
199                reduced_lin_reg.append(None)
200
201        sigs = [get_sig(r,orig_lin_reg,len(data)) for r in reduced_lin_reg]
202        if sigs and max(sigs) > remove_sig:
203            # remove that attribute, start again
204            crit_att = inc_atts[sigs.index(max(sigs))]
205            not_inc_atts.append(crit_att)
206            inc_atts.remove(crit_att)
207            changed_model = True
208            continue
209
210        # add all significant conditions (go through all attributes in not_inc_atts,
211        # is there one that significantly improves the model?
212        more_complex_lin_reg = []
213        for ati in range(len(not_inc_atts)):
214            try:
215                more_complex_lin_reg.append(LinearRegressionLearner(data, weight, use_attributes = inc_atts + [not_inc_atts[ati]]))
216            except:
217                more_complex_lin_reg.append(None)
218
219        sigs = [get_sig(orig_lin_reg,r,len(data)) for r in more_complex_lin_reg]
220        if sigs and min(sigs) < add_sig:
221            best_att = not_inc_atts[sigs.index(min(sigs))]
222            inc_atts.append(best_att)
223            not_inc_atts.remove(best_att)
224            changed_model = True
225    return inc_atts
226
227
228########################################################################
229# Partial Least-Squares Regression (PLS)
230
231# Function is used in PLSRegression
232def normalize(vector):
233    return vector / numpy.linalg.norm(vector)
234
235class PLSRegressionLearner(object):
236    """PLSRegressionLearner(data, y, x=None, nc=None)"""
237    def __new__(self, data=None, **kwds):
238        learner = object.__new__(self, **kwds)
239        if data:
240            learner.__init__(**kwds) # force init
241            return learner(data, **kwds)
242        else:
243            return learner  # invokes the __init__
244
245    def __init__(self, name='PLS regression', nc = None, save_partial=False, **kwds):
246        self.name = name
247        self.nc = nc
248        self.save_partial = save_partial #whether to save partial results
249
250    def __call__(self, data, y=None, x=None, nc=None, weight=None, **kwds):
251
252        if y == None:
253            y = [ data.domain.classVar ]
254        if x == None:
255            x = [v for v in data.domain.variables if v not in y]
256
257        Ncomp = nc if nc is not None else len(x)
258
259        dataX = orange.ExampleTable(orange.Domain(x, False), data)
260        dataY = orange.ExampleTable(orange.Domain(y, False), data)
261
262        # transformation to numpy arrays
263        X = dataX.toNumpy()[0]
264        Y = dataY.toNumpy()[0]
265
266        # data dimensions
267        n, mx = numpy.shape(X)
268        my = numpy.shape(Y)[1]
269
270        # Z-scores of original matrices
271        YMean = numpy.mean(Y, axis = 0)
272        YStd = numpy.std(Y, axis = 0)
273        XMean = numpy.mean(X, axis = 0)
274        XStd = numpy.std(X, axis = 0)
275
276        #FIXME: standard deviation should never be 0. Ask Lan, if the following
277        #fix is ok.
278        XStd = numpy.maximum(XStd, 10e-16)
279        YStd = numpy.maximum(YStd, 10e-16)
280
281        X = (X-XMean)/XStd
282        Y = (Y-YMean)/YStd
283
284        P = numpy.empty((mx,Ncomp))
285        C = numpy.empty((my,Ncomp))
286        T = numpy.empty((n,Ncomp))
287        U = numpy.empty((n,Ncomp))
288        B = numpy.zeros((Ncomp,Ncomp))
289        W = numpy.empty((mx,Ncomp))
290        E,F = X,Y
291
292        # main algorithm
293        for i in range(Ncomp):
294            u = numpy.random.random_sample((n,1)) #FIXME random seed?
295            w = normalize(dot(E.T,u))
296            t = normalize(dot(E,w))
297            dif = t
298
300            while numpy.linalg.norm(dif) > 10e-16:
301                c = normalize(dot(F.T,t))
302                u = dot(F,c)
303                w = normalize(dot(E.T,u))
304                t0 = normalize(dot(E,w))
305                dif = t - t0
306                t = t0
307
308            #print "T", T
309            #print "X*W", numpy.dot(X,W)
310
311            T[:,i] = t.T
312            U[:,i] = u.T
313            C[:,i] = c.T
314            W[:,i] = w.T
315            b = dot(t.T,u)[0]
316            B[i][i] = b
317            p = dot(E.T,t)
318            P[:,i] = p.T
319            E = E - dot(t,p.T)
320            F = F - b * dot(t,c.T)
321
322        # esimated Y
323        YE = dot(dot(T,B),C.T)*YStd + YMean
324        Y = Y*numpy.std(Y, axis = 0)+ YMean
325        BPls = dot(dot(numpy.linalg.pinv(P.T),B),C.T)
326
327        partial = {}
328        if self.save_partial:
329            partial["T"] = T
330            partial["U"] = U
331            partial["C"] = C
332            partial["W"] = W
333            partial["P"] = P
334
335        return PLSRegression(domain=data.domain, BPls=BPls, YMean=YMean, YStd=YStd, XMean=XMean, XStd=XStd, name=self.name, **partial)
336
337class PLSRegression(orange.Classifier):
338    def __init__(self, **kwargs):
339        for a,b in kwargs.items():
340            setattr(self, a, b)
341
342    def __call__(self, example):
343       ex = orange.Example(self.domain, example)
344       ex = numpy.array(ex.native()) #FIXME wrong -> no X and Y separation!
345       ex = (ex - self.XMean) / self.XStd
346       yhat = dot(ex, self.BPls) * self.YStd + self.YMean
347       return yhat
Note: See TracBrowser for help on using the repository browser.