source: orange/orange/Orange/classification/logreg.py @ 7747:b8cd550d6e36

Revision 7747:b8cd550d6e36, 42.2 KB checked in by jzbontar <jure.zbontar@…>, 3 years ago (diff)

Rename camel case identifiers in logistic regression.

Line 
1"""
2.. index: logistic regression
3.. index:
4   single: classification; logistic regression
5
6*******************
7Logistic regression
8*******************
9
10Implements `logistic regression
11<http://en.wikipedia.org/wiki/Logistic_regression>`_ with an extension for
12proper treatment of discrete features.  The algorithm can handle various
13anomalies in features, such as constant variables and singularities, that
14could make fitting of logistic regression almost impossible. Stepwise
15logistic regression, which iteratively selects the most informative
16features, is also supported.
17
18
19.. autoclass:: LogRegLearner
20.. autoclass:: StepWiseFSS
21.. autofunction:: dump
22
23Examples
24========
25
26The first example shows a very simple induction of a logistic regression
27classifier (`logreg-run.py`_, uses `titanic.tab`_).
28
29.. literalinclude:: code/logreg-run.py
30
31Result::
32
33    Classification accuracy: 0.778282598819
34
35    class attribute = survived
36    class values = <no, yes>
37
38        Attribute       beta  st. error     wald Z          P OR=exp(beta)
39
40        Intercept      -1.23       0.08     -15.15      -0.00
41     status=first       0.86       0.16       5.39       0.00       2.36
42    status=second      -0.16       0.18      -0.91       0.36       0.85
43     status=third      -0.92       0.15      -6.12       0.00       0.40
44        age=child       1.06       0.25       4.30       0.00       2.89
45       sex=female       2.42       0.14      17.04       0.00      11.25
46
47The next examples shows how to handle singularities in data sets
48(`logreg-singularities.py`_, uses `adult_sample.tab`_).
49
50.. literalinclude:: code/logreg-singularities.py
51
52The first few lines of the output of this script are::
53
54    <=50K <=50K
55    <=50K <=50K
56    <=50K <=50K
57    >50K >50K
58    <=50K >50K
59
60    class attribute = y
61    class values = <>50K, <=50K>
62
63                               Attribute       beta  st. error     wald Z          P OR=exp(beta)
64
65                               Intercept       6.62      -0.00       -inf       0.00
66                                     age      -0.04       0.00       -inf       0.00       0.96
67                                  fnlwgt      -0.00       0.00       -inf       0.00       1.00
68                           education-num      -0.28       0.00       -inf       0.00       0.76
69                 marital-status=Divorced       4.29       0.00        inf       0.00      72.62
70            marital-status=Never-married       3.79       0.00        inf       0.00      44.45
71                marital-status=Separated       3.46       0.00        inf       0.00      31.95
72                  marital-status=Widowed       3.85       0.00        inf       0.00      46.96
73    marital-status=Married-spouse-absent       3.98       0.00        inf       0.00      53.63
74        marital-status=Married-AF-spouse       4.01       0.00        inf       0.00      55.19
75                 occupation=Tech-support      -0.32       0.00       -inf       0.00       0.72
76
77If :obj:`removeSingular` is set to 0, inducing a logistic regression
78classifier would return an error::
79
80    Traceback (most recent call last):
81      File "logreg-singularities.py", line 4, in <module>
82        lr = classification.logreg.LogRegLearner(table, removeSingular=0)
83      File "/home/jure/devel/orange/Orange/classification/logreg.py", line 255, in LogRegLearner
84        return lr(examples, weightID)
85      File "/home/jure/devel/orange/Orange/classification/logreg.py", line 291, in __call__
86        lr = learner(examples, weight)
87    orange.KernelException: 'orange.LogRegLearner': singularity in workclass=Never-worked
88
89We can see that the attribute workclass is causing a singularity.
90
91The example below shows, how the use of stepwise logistic regression can help to
92gain in classification performance (`logreg-stepwise.py`_, uses `ionosphere.tab`_):
93
94.. literalinclude:: code/logreg-stepwise.py
95
96The output of this script is::
97
98    Learner      CA
99    logistic     0.841
100    filtered     0.846
101
102    Number of times attributes were used in cross-validation:
103     1 x a21
104    10 x a22
105     8 x a23
106     7 x a24
107     1 x a25
108    10 x a26
109    10 x a27
110     3 x a28
111     7 x a29
112     9 x a31
113     2 x a16
114     7 x a12
115     1 x a32
116     8 x a15
117    10 x a14
118     4 x a17
119     7 x a30
120    10 x a11
121     1 x a10
122     1 x a13
123    10 x a34
124     2 x a19
125     1 x a18
126    10 x a3
127    10 x a5
128     4 x a4
129     4 x a7
130     8 x a6
131    10 x a9
132    10 x a8
133
134.. _logreg-run.py: code/logreg-run.py
135.. _logreg-singularities.py: code/logreg-singularities.py
136.. _logreg-stepwise.py: code/logreg-stepwise.py
137
138.. _ionosphere.tab: code/ionosphere.tab
139.. _adult_sample.tab: code/adult_sample.tab
140.. _titanic.tab: code/titanic.tab
141
142"""
143
144#from Orange.core import LogRegLearner, LogRegClassifier, LogRegFitter, LogRegFitter_Cholesky
145
146import Orange
147import math, os
148import warnings
149from numpy import *
150from numpy.linalg import *
151
152
153##########################################################################
154## Print out methods
155
156def dump(classifier):
157    """ Formatted string of all major features in logistic
158    regression classifier.
159
160    :param classifier: logistic regression classifier
161    """
162
163    # print out class values
164    out = ['']
165    out.append("class attribute = " + classifier.domain.classVar.name)
166    out.append("class values = " + str(classifier.domain.classVar.values))
167    out.append('')
168   
169    # get the longest attribute name
170    longest=0
171    for at in classifier.continuizedDomain.attributes:
172        if len(at.name)>longest:
173            longest=len(at.name);
174
175    # print out the head
176    formatstr = "%"+str(longest)+"s %10s %10s %10s %10s %10s"
177    out.append(formatstr % ("Feature", "beta", "st. error", "wald Z", "P", "OR=exp(beta)"))
178    out.append('')
179    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f"   
180    out.append(formatstr % ("Intercept", classifier.beta[0], classifier.beta_se[0], classifier.wald_Z[0], classifier.P[0]))
181    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f %10.2f"   
182    for i in range(len(classifier.continuizedDomain.attributes)):
183        out.append(formatstr % (classifier.continuizedDomain.attributes[i].name, classifier.beta[i+1], classifier.beta_se[i+1], classifier.wald_Z[i+1], abs(classifier.P[i+1]), math.exp(classifier.beta[i+1])))
184
185    return '\n'.join(out)
186       
187
188def has_discrete_values(domain):
189    for at in domain.attributes:
190        if at.varType == Orange.core.VarTypes.Discrete:
191            return 1
192    return 0
193
194class LogRegLearner(Orange.classification.Learner):
195    """ Logistic regression learner.
196
197    Implements logistic regression. If data instances are provided to
198    the constructor, the learning algorithm is called and the resulting
199    classifier is returned instead of the learner.
200
201    :param table: data table with either discrete or continuous features
202    :type table: Orange.data.Table
203    :param weightID: the ID of the weight meta attribute
204    :type weightID: int
205    :param removeSingular: set to 1 if you want automatic removal of disturbing features, such as constants and singularities
206    :type removeSingular: bool
207    :param fitter: the fitting algorithm (by default the Newton-Raphson fitting algorithm is used)
208    :param stepwiseLR: set to 1 if you wish to use stepwise logistic regression
209    :type stepwiseLR: bool
210    :param addCrit: parameter for stepwise feature selection
211    :type addCrit: float
212    :param deleteCrit: parameter for stepwise feature selection
213    :type deleteCrit: float
214    :param numFeatures: parameter for stepwise feature selection
215    :type numFeatures: int
216    :rtype: :obj:`LogRegLearner` or :obj:`LogRegClassifier`
217
218    """
219    def __new__(cls, instances=None, weightID=0, **argkw):
220        self = Orange.classification.Learner.__new__(cls, **argkw)
221        if instances:
222            self.__init__(**argkw)
223            return self.__call__(instances, weightID)
224        else:
225            return self
226
227    def __init__(self, removeSingular=0, fitter = None, **kwds):
228        self.__dict__.update(kwds)
229        self.removeSingular = removeSingular
230        self.fitter = None
231
232    def __call__(self, examples, weight=0):
233        imputer = getattr(self, "imputer", None) or None
234        if getattr(self, "removeMissing", 0):
235            examples = Orange.core.Preprocessor_dropMissing(examples)
236##        if hasDiscreteValues(examples.domain):
237##            examples = createNoDiscTable(examples)
238        if not len(examples):
239            return None
240        if getattr(self, "stepwiseLR", 0):
241            addCrit = getattr(self, "addCrit", 0.2)
242            removeCrit = getattr(self, "removeCrit", 0.3)
243            numFeatures = getattr(self, "numFeatures", -1)
244            attributes = StepWiseFSS(examples, addCrit = addCrit, deleteCrit = removeCrit, imputer = imputer, numFeatures = numFeatures)
245            tmpDomain = Orange.core.Domain(attributes, examples.domain.classVar)
246            tmpDomain.addmetas(examples.domain.getmetas())
247            examples = examples.select(tmpDomain)
248        learner = Orange.core.LogRegLearner()
249        learner.imputerConstructor = imputer
250        if imputer:
251            examples = self.imputer(examples)(examples)
252        examples = Orange.core.Preprocessor_dropMissing(examples)
253        if self.fitter:
254            learner.fitter = self.fitter
255        if self.removeSingular:
256            lr = learner.fitModel(examples, weight)
257        else:
258            lr = learner(examples, weight)
259        while isinstance(lr, Orange.core.Variable):
260            if isinstance(lr.getValueFrom, Orange.core.ClassifierFromVar) and isinstance(lr.getValueFrom.transformer, Orange.core.Discrete2Continuous):
261                lr = lr.getValueFrom.variable
262            attributes = examples.domain.attributes[:]
263            if lr in attributes:
264                attributes.remove(lr)
265            else:
266                attributes.remove(lr.getValueFrom.variable)
267            newDomain = Orange.core.Domain(attributes, examples.domain.classVar)
268            newDomain.addmetas(examples.domain.getmetas())
269            examples = examples.select(newDomain)
270            lr = learner.fitModel(examples, weight)
271        return lr
272
273
274
275class UnivariateLogRegLearner(Orange.classification.Learner):
276    def __new__(cls, instances=None, **argkw):
277        self = Orange.classification.Learner.__new__(cls, **argkw)
278        if instances:
279            self.__init__(**argkw)
280            return self.__call__(instances)
281        else:
282            return self
283
284    def __init__(self, **kwds):
285        self.__dict__.update(kwds)
286
287    def __call__(self, examples):
288        examples = createFullNoDiscTable(examples)
289        classifiers = map(lambda x: LogRegLearner(Orange.core.Preprocessor_dropMissing(examples.select(Orange.core.Domain(x, examples.domain.classVar)))), examples.domain.attributes)
290        maj_classifier = LogRegLearner(Orange.core.Preprocessor_dropMissing(examples.select(Orange.core.Domain(examples.domain.classVar))))
291        beta = [maj_classifier.beta[0]] + [x.beta[1] for x in classifiers]
292        beta_se = [maj_classifier.beta_se[0]] + [x.beta_se[1] for x in classifiers]
293        P = [maj_classifier.P[0]] + [x.P[1] for x in classifiers]
294        wald_Z = [maj_classifier.wald_Z[0]] + [x.wald_Z[1] for x in classifiers]
295        domain = examples.domain
296
297        return Univariate_LogRegClassifier(beta = beta, beta_se = beta_se, P = P, wald_Z = wald_Z, domain = domain)
298
299class UnivariateLogRegClassifier(Orange.core.Classifier):
300    def __init__(self, **kwds):
301        self.__dict__.update(kwds)
302
303    def __call__(self, example, resultType = Orange.core.GetValue):
304        # classification not implemented yet. For now its use is only to provide regression coefficients and its statistics
305        pass
306   
307
308class LogRegLearnerGetPriors(object):
309    def __new__(cls, instances=None, weightID=0, **argkw):
310        self = object.__new__(cls, **argkw)
311        if instances:
312            self.__init__(**argkw)
313            return self.__call__(instances, weightID)
314        else:
315            return self
316
317    def __init__(self, removeSingular=0, **kwds):
318        self.__dict__.update(kwds)
319        self.removeSingular = removeSingular
320    def __call__(self, examples, weight=0):
321        # next function changes data set to a extended with unknown values
322        def createLogRegExampleTable(data, weightID):
323            setsOfData = []
324            for at in data.domain.attributes:
325                # za vsak atribut kreiraj nov newExampleTable newData
326                # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable
327                if at.varType == Orange.core.VarTypes.Continuous:
328                    atDisc = Orange.core.FloatVariable(at.name + "Disc")
329                    newDomain = Orange.core.Domain(data.domain.attributes+[atDisc,data.domain.classVar])
330                    newDomain.addmetas(data.domain.getmetas())
331                    newData = Orange.core.ExampleTable(newDomain,data)
332                    altData = Orange.core.ExampleTable(newDomain,data)
333                    for i,d in enumerate(newData):
334                        d[atDisc] = 0
335                        d[weightID] = 1*data[i][weightID]
336                    for i,d in enumerate(altData):
337                        d[atDisc] = 1
338                        d[at] = 0
339                        d[weightID] = 0.000001*data[i][weightID]
340                elif at.varType == Orange.core.VarTypes.Discrete:
341                # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X"
342                    atNew = Orange.core.EnumVariable(at.name, values = at.values + [at.name+"X"])
343                    newDomain = Orange.core.Domain(filter(lambda x: x!=at, data.domain.attributes)+[atNew,data.domain.classVar])
344                    newDomain.addmetas(data.domain.getmetas())
345                    newData = Orange.core.ExampleTable(newDomain,data)
346                    altData = Orange.core.ExampleTable(newDomain,data)
347                    for i,d in enumerate(newData):
348                        d[atNew] = data[i][at]
349                        d[weightID] = 1*data[i][weightID]
350                    for i,d in enumerate(altData):
351                        d[atNew] = at.name+"X"
352                        d[weightID] = 0.000001*data[i][weightID]
353                newData.extend(altData)
354                setsOfData.append(newData)
355            return setsOfData
356                 
357        learner = LogRegLearner(imputer = Orange.core.ImputerConstructor_average(), removeSingular = self.removeSingular)
358        # get Original Model
359        orig_model = learner(examples,weight)
360        if orig_model.fit_status:
361            print "Warning: model did not converge"
362
363        # get extended Model (you should not change data)
364        if weight == 0:
365            weight = Orange.core.newmetaid()
366            examples.addMetaAttribute(weight, 1.0)
367        extended_set_of_examples = createLogRegExampleTable(examples, weight)
368        extended_models = [learner(extended_examples, weight) \
369                           for extended_examples in extended_set_of_examples]
370
371##        print examples[0]
372##        printOUT(orig_model)
373##        print orig_model.domain
374##        print orig_model.beta
375##        print orig_model.beta[orig_model.continuizedDomain.attributes[-1]]
376##        for i,m in enumerate(extended_models):
377##            print examples.domain.attributes[i]
378##            printOUT(m)
379           
380       
381        # izracunas odstopanja
382        # get sum of all betas
383        beta = 0
384        betas_ap = []
385        for m in extended_models:
386            beta_add = m.beta[m.continuizedDomain.attributes[-1]]
387            betas_ap.append(beta_add)
388            beta = beta + beta_add
389       
390        # substract it from intercept
391        #print "beta", beta
392        logistic_prior = orig_model.beta[0]+beta
393       
394        # compare it to bayes prior
395        bayes = Orange.core.BayesLearner(examples)
396        bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0])
397
398        # normalize errors
399##        print "bayes", bayes_prior
400##        print "lr", orig_model.beta[0]
401##        print "lr2", logistic_prior
402##        print "dist", Orange.core.Distribution(examples.domain.classVar,examples)
403##        print "prej", betas_ap
404
405        # error normalization - to avoid errors due to assumption of independence of unknown values
406        dif = bayes_prior - logistic_prior
407        positives = sum(filter(lambda x: x>=0, betas_ap))
408        negatives = -sum(filter(lambda x: x<0, betas_ap))
409        if not negatives == 0:
410            kPN = positives/negatives
411            diffNegatives = dif/(1+kPN)
412            diffPositives = kPN*diffNegatives
413            kNegatives = (negatives-diffNegatives)/negatives
414            kPositives = positives/(positives-diffPositives)
415    ##        print kNegatives
416    ##        print kPositives
417
418            for i,b in enumerate(betas_ap):
419                if b<0: betas_ap[i]*=kNegatives
420                else: betas_ap[i]*=kPositives
421        #print "potem", betas_ap
422
423        # vrni originalni model in pripadajoce apriorne niclele
424        return (orig_model, betas_ap)
425        #return (bayes_prior,orig_model.beta[examples.domain.classVar],logistic_prior)
426
427class LogRegLearnerGetPriorsOneTable:
428    def __init__(self, removeSingular=0, **kwds):
429        self.__dict__.update(kwds)
430        self.removeSingular = removeSingular
431    def __call__(self, examples, weight=0):
432        # next function changes data set to a extended with unknown values
433        def createLogRegExampleTable(data, weightID):
434            finalData = Orange.core.ExampleTable(data)
435            origData = Orange.core.ExampleTable(data)
436            for at in data.domain.attributes:
437                # za vsak atribut kreiraj nov newExampleTable newData
438                # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable
439                if at.varType == Orange.core.VarTypes.Continuous:
440                    atDisc = Orange.core.FloatVariable(at.name + "Disc")
441                    newDomain = Orange.core.Domain(origData.domain.attributes+[atDisc,data.domain.classVar])
442                    newDomain.addmetas(newData.domain.getmetas())
443                    finalData = Orange.core.ExampleTable(newDomain,finalData)
444                    newData = Orange.core.ExampleTable(newDomain,origData)
445                    origData = Orange.core.ExampleTable(newDomain,origData)
446                    for d in origData:
447                        d[atDisc] = 0
448                    for d in finalData:
449                        d[atDisc] = 0
450                    for i,d in enumerate(newData):
451                        d[atDisc] = 1
452                        d[at] = 0
453                        d[weightID] = 100*data[i][weightID]
454                       
455                elif at.varType == Orange.core.VarTypes.Discrete:
456                # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X"
457                    atNew = Orange.core.EnumVariable(at.name, values = at.values + [at.name+"X"])
458                    newDomain = Orange.core.Domain(filter(lambda x: x!=at, origData.domain.attributes)+[atNew,origData.domain.classVar])
459                    newDomain.addmetas(origData.domain.getmetas())
460                    temp_finalData = Orange.core.ExampleTable(finalData)
461                    finalData = Orange.core.ExampleTable(newDomain,finalData)
462                    newData = Orange.core.ExampleTable(newDomain,origData)
463                    temp_origData = Orange.core.ExampleTable(origData)
464                    origData = Orange.core.ExampleTable(newDomain,origData)
465                    for i,d in enumerate(origData):
466                        d[atNew] = temp_origData[i][at]
467                    for i,d in enumerate(finalData):
468                        d[atNew] = temp_finalData[i][at]                       
469                    for i,d in enumerate(newData):
470                        d[atNew] = at.name+"X"
471                        d[weightID] = 10*data[i][weightID]
472                finalData.extend(newData)
473            return finalData
474                 
475        learner = LogRegLearner(imputer = Orange.core.ImputerConstructor_average(), removeSingular = self.removeSingular)
476        # get Original Model
477        orig_model = learner(examples,weight)
478
479        # get extended Model (you should not change data)
480        if weight == 0:
481            weight = Orange.core.newmetaid()
482            examples.addMetaAttribute(weight, 1.0)
483        extended_examples = createLogRegExampleTable(examples, weight)
484        extended_model = learner(extended_examples, weight)
485
486##        print examples[0]
487##        printOUT(orig_model)
488##        print orig_model.domain
489##        print orig_model.beta
490
491##        printOUT(extended_model)       
492        # izracunas odstopanja
493        # get sum of all betas
494        beta = 0
495        betas_ap = []
496        for m in extended_models:
497            beta_add = m.beta[m.continuizedDomain.attributes[-1]]
498            betas_ap.append(beta_add)
499            beta = beta + beta_add
500       
501        # substract it from intercept
502        #print "beta", beta
503        logistic_prior = orig_model.beta[0]+beta
504       
505        # compare it to bayes prior
506        bayes = Orange.core.BayesLearner(examples)
507        bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0])
508
509        # normalize errors
510        #print "bayes", bayes_prior
511        #print "lr", orig_model.beta[0]
512        #print "lr2", logistic_prior
513        #print "dist", Orange.core.Distribution(examples.domain.classVar,examples)
514        k = (bayes_prior-orig_model.beta[0])/(logistic_prior-orig_model.beta[0])
515        #print "prej", betas_ap
516        betas_ap = [k*x for x in betas_ap]               
517        #print "potem", betas_ap
518
519        # vrni originalni model in pripadajoce apriorne niclele
520        return (orig_model, betas_ap)
521        #return (bayes_prior,orig_model.beta[data.domain.classVar],logistic_prior)
522
523
524######################################
525#### Fitters for logistic regression (logreg) learner ####
526######################################
527
528def pr(x, betas):
529    k = math.exp(dot(x, betas))
530    return k / (1+k)
531
532def lh(x,y,betas):
533    llh = 0.0
534    for i,x_i in enumerate(x):
535        pr = pr(x_i,betas)
536        llh += y[i]*log(max(pr,1e-6)) + (1-y[i])*log(max(1-pr,1e-6))
537    return llh
538
539
540def diag(vector):
541    mat = identity(len(vector), Float)
542    for i,v in enumerate(vector):
543        mat[i][i] = v
544    return mat
545   
546class SimpleFitter(Orange.core.LogRegFitter):
547    def __init__(self, penalty=0, se_penalty = False):
548        self.penalty = penalty
549        self.se_penalty = se_penalty
550    def __call__(self, data, weight=0):
551        ml = data.native(0)
552        for i in range(len(data.domain.attributes)):
553          a = data.domain.attributes[i]
554          if a.varType == Orange.core.VarTypes.Discrete:
555            for m in ml:
556              m[i] = a.values.index(m[i])
557        for m in ml:
558          m[-1] = data.domain.classVar.values.index(m[-1])
559        Xtmp = array(ml)
560        y = Xtmp[:,-1]   # true probabilities (1's or 0's)
561        one = reshape(array([1]*len(data)), (len(data),1)) # intercept column
562        X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data
563
564        betas = array([0.0] * (len(data.domain.attributes)+1))
565        oldBetas = array([1.0] * (len(data.domain.attributes)+1))
566        N = len(data)
567
568        pen_matrix = array([self.penalty] * (len(data.domain.attributes)+1))
569        if self.se_penalty:
570            p = array([pr(X[i], betas) for i in range(len(data))])
571            W = identity(len(data), Float)
572            pp = p * (1.0-p)
573            for i in range(N):
574                W[i,i] = pp[i]
575            se = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X)))))
576            for i,p in enumerate(pen_matrix):
577                pen_matrix[i] *= se[i]
578        # predict the probability for an instance, x and betas are vectors
579        # start the computation
580        likelihood = 0.
581        likelihood_new = 1.
582        while abs(likelihood - likelihood_new)>1e-5:
583            likelihood = likelihood_new
584            oldBetas = betas
585            p = array([pr(X[i], betas) for i in range(len(data))])
586
587            W = identity(len(data), Float)
588            pp = p * (1.0-p)
589            for i in range(N):
590                W[i,i] = pp[i]
591
592            WI = inverse(W)
593            z = matrixmultiply(X, betas) + matrixmultiply(WI, y - p)
594
595            tmpA = inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))+diag(pen_matrix))
596            tmpB = matrixmultiply(transpose(X), y-p)
597            betas = oldBetas + matrixmultiply(tmpA,tmpB)
598#            betaTemp = matrixmultiply(matrixmultiply(matrixmultiply(matrixmultiply(tmpA,transpose(X)),W),X),oldBetas)
599#            print betaTemp
600#            tmpB = matrixmultiply(transpose(X), matrixmultiply(W, z))
601#            betas = matrixmultiply(tmpA, tmpB)
602            likelihood_new = lh(X,y,betas)-self.penalty*sum([b*b for b in betas])
603            print likelihood_new
604
605           
606           
607##        XX = sqrt(diagonal(inverse(matrixmultiply(transpose(X),X))))
608##        yhat = array([pr(X[i], betas) for i in range(len(data))])
609##        ss = sum((y - yhat) ** 2) / (N - len(data.domain.attributes) - 1)
610##        sigma = math.sqrt(ss)
611        p = array([pr(X[i], betas) for i in range(len(data))])
612        W = identity(len(data), Float)
613        pp = p * (1.0-p)
614        for i in range(N):
615            W[i,i] = pp[i]
616        diXWX = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X)))))
617        xTemp = matrixmultiply(matrixmultiply(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))),transpose(X)),y)
618        beta = []
619        beta_se = []
620        print "likelihood ridge", likelihood
621        for i in range(len(betas)):
622            beta.append(betas[i])
623            beta_se.append(diXWX[i])
624        return (self.OK, beta, beta_se, 0)
625
626def pr_bx(bx):
627    if bx > 35:
628        return 1
629    if bx < -35:
630        return 0
631    return exp(bx)/(1+exp(bx))
632
633class BayesianFitter(Orange.core.LogRegFitter):
634    def __init__(self, penalty=0, anch_examples=[], tau = 0):
635        self.penalty = penalty
636        self.anch_examples = anch_examples
637        self.tau = tau
638
639    def create_array_data(self,data):
640        if not len(data):
641            return (array([]),array([]))
642        # convert data to numeric
643        ml = data.native(0)
644        for i,a in enumerate(data.domain.attributes):
645          if a.varType == Orange.core.VarTypes.Discrete:
646            for m in ml:
647              m[i] = a.values.index(m[i])
648        for m in ml:
649          m[-1] = data.domain.classVar.values.index(m[-1])
650        Xtmp = array(ml)
651        y = Xtmp[:,-1]   # true probabilities (1's or 0's)
652        one = reshape(array([1]*len(data)), (len(data),1)) # intercept column
653        X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data
654        return (X,y)
655   
656    def __call__(self, data, weight=0):
657        (X,y)=self.create_array_data(data)
658
659        exTable = Orange.core.ExampleTable(data.domain)
660        for id,ex in self.anch_examples:
661            exTable.extend(Orange.core.ExampleTable(ex,data.domain))
662        (X_anch,y_anch)=self.create_array_data(exTable)
663
664        betas = array([0.0] * (len(data.domain.attributes)+1))
665
666        likelihood,betas = self.estimate_beta(X,y,betas,[0]*(len(betas)),X_anch,y_anch)
667
668        # get attribute groups atGroup = [(startIndex, number of values), ...)
669        ats = data.domain.attributes
670        atVec=reduce(lambda x,y: x+[(y,not y==x[-1][0])], [a.getValueFrom and a.getValueFrom.whichVar or a for a in ats],[(ats[0].getValueFrom and ats[0].getValueFrom.whichVar or ats[0],0)])[1:]
671        atGroup=[[0,0]]
672        for v_i,v in enumerate(atVec):
673            if v[1]==0: atGroup[-1][1]+=1
674            else:       atGroup.append([v_i,1])
675       
676        # compute zero values for attributes
677        sumB = 0.
678        for ag in atGroup:
679            X_temp = concatenate((X[:,:ag[0]+1],X[:,ag[0]+1+ag[1]:]),1)
680            if X_anch:
681                X_anch_temp = concatenate((X_anch[:,:ag[0]+1],X_anch[:,ag[0]+1+ag[1]:]),1)
682            else: X_anch_temp = X_anch
683##            print "1", concatenate((betas[:i+1],betas[i+2:]))
684##            print "2", betas
685            likelihood_temp,betas_temp=self.estimate_beta(X_temp,y,concatenate((betas[:ag[0]+1],betas[ag[0]+ag[1]+1:])),[0]+[1]*(len(betas)-1-ag[1]),X_anch_temp,y_anch)
686            print "finBetas", betas, betas_temp
687            print "betas", betas[0], betas_temp[0]
688            sumB += betas[0]-betas_temp[0]
689        apriori = Orange.core.Distribution(data.domain.classVar, data)
690        aprioriProb = apriori[0]/apriori.abs
691       
692        print "koncni rezultat", sumB, math.log((1-aprioriProb)/aprioriProb), betas[0]
693           
694        beta = []
695        beta_se = []
696        print "likelihood2", likelihood
697        for i in range(len(betas)):
698            beta.append(betas[i])
699            beta_se.append(0.0)
700        return (self.OK, beta, beta_se, 0)
701
702     
703       
704    def estimate_beta(self,X,y,betas,const_betas,X_anch,y_anch):
705        N,N_anch = len(y),len(y_anch)
706        r,r_anch = array([dot(X[i], betas) for i in range(N)]),\
707                   array([dot(X_anch[i], betas) for i in range(N_anch)])
708        p    = array([pr_bx(ri) for ri in r])
709        X_sq = X*X
710
711        max_delta      = [1.]*len(const_betas)
712        likelihood     = -1.e+10
713        likelihood_new = -1.e+9
714        while abs(likelihood - likelihood_new)>0.01 and max(max_delta)>0.01:
715            likelihood = likelihood_new
716            print likelihood
717            betas_temp = [b for b in betas]
718            for j in range(len(betas)):
719                if const_betas[j]: continue
720                dl = matrixmultiply(X[:,j],transpose(y-p))
721                for xi,x in enumerate(X_anch):
722                    dl += self.penalty*x[j]*(y_anch[xi] - pr_bx(r_anch[xi]*self.penalty))
723
724                ddl = matrixmultiply(X_sq[:,j],transpose(p*(1-p)))
725                for xi,x in enumerate(X_anch):
726                    ddl += self.penalty*x[j]*pr_bx(r[xi]*self.penalty)*(1-pr_bx(r[xi]*self.penalty))
727
728                if j==0:
729                    dv = dl/max(ddl,1e-6)
730                elif betas[j] == 0: # special handling due to non-defined first and second derivatives
731                    dv = (dl-self.tau)/max(ddl,1e-6)
732                    if dv < 0:
733                        dv = (dl+self.tau)/max(ddl,1e-6)
734                        if dv > 0:
735                            dv = 0
736                else:
737                    dl -= sign(betas[j])*self.tau
738                    dv = dl/max(ddl,1e-6)
739                    if not sign(betas[j] + dv) == sign(betas[j]):
740                        dv = -betas[j]
741                dv = min(max(dv,-max_delta[j]),max_delta[j])
742                r+= X[:,j]*dv
743                p = array([pr_bx(ri) for ri in r])
744                if N_anch:
745                    r_anch+=X_anch[:,j]*dv
746                betas[j] += dv
747                max_delta[j] = max(2*abs(dv),max_delta[j]/2)
748            likelihood_new = lh(X,y,betas)
749            for xi,x in enumerate(X_anch):
750                try:
751                    likelihood_new += y_anch[xi]*r_anch[xi]*self.penalty-log(1+exp(r_anch[xi]*self.penalty))
752                except:
753                    likelihood_new += r_anch[xi]*self.penalty*(y_anch[xi]-1)
754            likelihood_new -= sum([abs(b) for b in betas[1:]])*self.tau
755            if likelihood_new < likelihood:
756                max_delta = [md/4 for md in max_delta]
757                likelihood_new = likelihood
758                likelihood = likelihood_new + 1.
759                betas = [b for b in betas_temp]
760        print "betas", betas
761        print "init_like", likelihood_new
762        print "pure_like", lh(X,y,betas)
763        return (likelihood,betas)
764   
765############################################################
766#  Feature subset selection for logistic regression
767
768def get_likelihood(fitter, examples):
769    res = fitter(examples)
770    if res[0] in [fitter.OK]: #, fitter.Infinity, fitter.Divergence]:
771       status, beta, beta_se, likelihood = res
772       if sum([abs(b) for b in beta])<sum([abs(b) for b in beta_se]):
773           return -100*len(examples)
774       return likelihood
775    else:
776       return -100*len(examples)
777       
778
779
780class StepWiseFSS(Orange.classification.Learner):
781  """Implementation of algorithm described in [Hosmer and Lemeshow, Applied Logistic Regression, 2000].
782
783  Perform stepwise logistic regression and return a list of the
784  most "informative" features. Each step of the algorithm is composed
785  of two parts. The first is backward elimination, where each already
786  chosen feature is tested for a significant contribution to the overall
787  model. If the worst among all tested features has higher significance
788  than is specified in :obj:`deleteCrit`, the feature is removed from
789  the model. The second step is forward selection, which is similar to
790  backward elimination. It loops through all the features that are not
791  in the model and tests whether they contribute to the common model
792  with significance lower that :obj:`addCrit`. The algorithm stops when
793  no feature in the model is to be removed and no feature not in the
794  model is to be added. By setting :obj:`numFeatures` larger than -1,
795  the algorithm will stop its execution when the number of features in model
796  exceeds that number.
797
798  Significances are assesed via the likelihood ration chi-square
799  test. Normal F test is not appropriate, because errors are assumed to
800  follow a binomial distribution.
801
802  If :obj:`table` is specified, stepwise logistic regression implemented
803  in :obj:`StepWiseFSS` is performed and a list of chosen features
804  is returned. If :obj:`table` is not specified an instance of
805  :obj:`StepWiseFSS` with all parameters set is returned.
806
807  :param table: data set
808  :type table: Orange.data.Table
809
810  :param addCrit: "Alpha" level to judge if variable has enough importance to be added in the new set. (e.g. if addCrit is 0.2, then features is added if its P is lower than 0.2)
811  :type addCrit: float
812
813  :param deleteCrit: Similar to addCrit, just that it is used at backward elimination. It should be higher than addCrit!
814  :type deleteCrit: float
815
816  :param numFeatures: maximum number of selected features, use -1 for infinity.
817  :type numFeatures: int
818  :rtype: :obj:`StepWiseFSS` or list of features
819
820  """
821
822  def __new__(cls, instances=None, **argkw):
823      self = Orange.classification.Learner.__new__(cls, **argkw)
824      if instances:
825          self.__init__(**argkw)
826          return self.__call__(instances)
827      else:
828          return self
829
830
831  def __init__(self, addCrit=0.2, deleteCrit=0.3, numFeatures = -1, **kwds):
832    self.__dict__.update(kwds)
833    self.addCrit = addCrit
834    self.deleteCrit = deleteCrit
835    self.numFeatures = numFeatures
836  def __call__(self, examples):
837    if getattr(self, "imputer", 0):
838        examples = self.imputer(examples)(examples)
839    if getattr(self, "removeMissing", 0):
840        examples = Orange.core.Preprocessor_dropMissing(examples)
841    continuizer = Orange.core.DomainContinuizer(zeroBased=1,continuousTreatment=Orange.core.DomainContinuizer.Leave,
842                                           multinomialTreatment = Orange.core.DomainContinuizer.FrequentIsBase,
843                                           classTreatment = Orange.core.DomainContinuizer.Ignore)
844    attr = []
845    remain_attr = examples.domain.attributes[:]
846
847    # get LL for Majority Learner
848    tempDomain = Orange.core.Domain(attr,examples.domain.classVar)
849    #tempData  = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
850    tempData  = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
851
852    ll_Old = get_likelihood(Orange.core.LogRegFitter_Cholesky(), tempData)
853    ll_Best = -1000000
854    length_Old = float(len(tempData))
855
856    stop = 0
857    while not stop:
858        # LOOP until all variables are added or no further deletion nor addition of attribute is possible
859        worstAt = None
860        # if there are more than 1 attribute then perform backward elimination
861        if len(attr) >= 2:
862            minG = 1000
863            worstAt = attr[0]
864            ll_Best = ll_Old
865            length_Best = length_Old
866            for at in attr:
867                # check all attribute whether its presence enough increases LL?
868
869                tempAttr = filter(lambda x: x!=at, attr)
870                tempDomain = Orange.core.Domain(tempAttr,examples.domain.classVar)
871                tempDomain.addmetas(examples.domain.getmetas())
872                # domain, calculate P for LL improvement.
873                tempDomain  = continuizer(Orange.core.Preprocessor_dropMissing(examples.select(tempDomain)))
874                tempData = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
875
876                ll_Delete = get_likelihood(Orange.core.LogRegFitter_Cholesky(), tempData)
877                length_Delete = float(len(tempData))
878                length_Avg = (length_Delete + length_Old)/2.0
879
880                G=-2*length_Avg*(ll_Delete/length_Delete-ll_Old/length_Old)
881
882                # set new worst attribute               
883                if G<minG:
884                    worstAt = at
885                    minG=G
886                    ll_Best = ll_Delete
887                    length_Best = length_Delete
888            # deletion of attribute
889           
890            if worstAt.varType==Orange.core.VarTypes.Continuous:
891                P=lchisqprob(minG,1);
892            else:
893                P=lchisqprob(minG,len(worstAt.values)-1);
894            if P>=self.deleteCrit:
895                attr.remove(worstAt)
896                remain_attr.append(worstAt)
897                nodeletion=0
898                ll_Old = ll_Best
899                length_Old = length_Best
900            else:
901                nodeletion=1
902        else:
903            nodeletion = 1
904            # END OF DELETION PART
905           
906        # if enough attributes has been chosen, stop the procedure
907        if self.numFeatures>-1 and len(attr)>=self.numFeatures:
908            remain_attr=[]
909         
910        # for each attribute in the remaining
911        maxG=-1
912        ll_Best = ll_Old
913        length_Best = length_Old
914        bestAt = None
915        for at in remain_attr:
916            tempAttr = attr + [at]
917            tempDomain = Orange.core.Domain(tempAttr,examples.domain.classVar)
918            tempDomain.addmetas(examples.domain.getmetas())
919            # domain, calculate P for LL improvement.
920            tempDomain  = continuizer(Orange.core.Preprocessor_dropMissing(examples.select(tempDomain)))
921            tempData = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
922            ll_New = get_likelihood(Orange.core.LogRegFitter_Cholesky(), tempData)
923
924            length_New = float(len(tempData)) # get number of examples in tempData to normalize likelihood
925
926            # P=PR(CHI^2>G), G=-2(L(0)-L(1))=2(E(0)-E(1))
927            length_avg = (length_New + length_Old)/2
928            G=-2*length_avg*(ll_Old/length_Old-ll_New/length_New);
929            if G>maxG:
930                bestAt = at
931                maxG=G
932                ll_Best = ll_New
933                length_Best = length_New
934        if not bestAt:
935            stop = 1
936            continue
937       
938        if bestAt.varType==Orange.core.VarTypes.Continuous:
939            P=lchisqprob(maxG,1);
940        else:
941            P=lchisqprob(maxG,len(bestAt.values)-1);
942        # Add attribute with smallest P to attributes(attr)
943        if P<=self.addCrit:
944            attr.append(bestAt)
945            remain_attr.remove(bestAt)
946            ll_Old = ll_Best
947            length_Old = length_Best
948
949        if (P>self.addCrit and nodeletion) or (bestAt == worstAt):
950            stop = 1
951
952    return attr
953
954
955class StepWiseFSSFilter(object):
956    def __new__(cls, instances=None, **argkw):
957        self = object.__new__(cls, **argkw)
958        if instances:
959            self.__init__(**argkw)
960            return self.__call__(instances)
961        else:
962            return self
963   
964    def __init__(self, addCrit=0.2, deleteCrit=0.3, numFeatures = -1):
965        self.addCrit = addCrit
966        self.deleteCrit = deleteCrit
967        self.numFeatures = numFeatures
968
969    def __call__(self, examples):
970        attr = StepWiseFSS(examples, addCrit=self.addCrit, deleteCrit = self.deleteCrit, numFeatures = self.numFeatures)
971        return examples.select(Orange.core.Domain(attr, examples.domain.classVar))
972               
973
974####################################
975##  PROBABILITY CALCULATIONS
976
977def lchisqprob(chisq,df):
978    """
979    Return the (1-tailed) probability value associated with the provided
980    chi-square value and df.  Adapted from chisq.c in Gary Perlman's |Stat.
981    """
982    BIG = 20.0
983    def ex(x):
984        BIG = 20.0
985        if x < -BIG:
986            return 0.0
987        else:
988            return math.exp(x)
989    if chisq <=0 or df < 1:
990        return 1.0
991    a = 0.5 * chisq
992    if df%2 == 0:
993        even = 1
994    else:
995        even = 0
996    if df > 1:
997        y = ex(-a)
998    if even:
999        s = y
1000    else:
1001        s = 2.0 * zprob(-math.sqrt(chisq))
1002    if (df > 2):
1003        chisq = 0.5 * (df - 1.0)
1004        if even:
1005            z = 1.0
1006        else:
1007            z = 0.5
1008        if a > BIG:
1009            if even:
1010                e = 0.0
1011            else:
1012                e = math.log(math.sqrt(math.pi))
1013            c = math.log(a)
1014            while (z <= chisq):
1015                e = math.log(z) + e
1016                s = s + ex(c*z-a-e)
1017                z = z + 1.0
1018            return s
1019        else:
1020            if even:
1021                e = 1.0
1022            else:
1023                e = 1.0 / math.sqrt(math.pi) / math.sqrt(a)
1024            c = 0.0
1025            while (z <= chisq):
1026                e = e * (a/float(z))
1027                c = c + e
1028                z = z + 1.0
1029            return (c*y+s)
1030    else:
1031        return s
1032
1033
1034def zprob(z):
1035    """
1036    Returns the area under the normal curve 'to the left of' the given z value.
1037    Thus::
1038
1039    for z<0, zprob(z) = 1-tail probability
1040    for z>0, 1.0-zprob(z) = 1-tail probability
1041    for any z, 2.0*(1.0-zprob(abs(z))) = 2-tail probability
1042
1043    Adapted from z.c in Gary Perlman's |Stat.
1044    """
1045    Z_MAX = 6.0    # maximum meaningful z-value
1046    if z == 0.0:
1047    x = 0.0
1048    else:
1049    y = 0.5 * math.fabs(z)
1050    if y >= (Z_MAX*0.5):
1051        x = 1.0
1052    elif (y < 1.0):
1053        w = y*y
1054        x = ((((((((0.000124818987 * w
1055            -0.001075204047) * w +0.005198775019) * w
1056              -0.019198292004) * w +0.059054035642) * w
1057            -0.151968751364) * w +0.319152932694) * w
1058          -0.531923007300) * w +0.797884560593) * y * 2.0
1059    else:
1060        y = y - 2.0
1061        x = (((((((((((((-0.000045255659 * y
1062                 +0.000152529290) * y -0.000019538132) * y
1063               -0.000676904986) * y +0.001390604284) * y
1064             -0.000794620820) * y -0.002034254874) * y
1065               +0.006549791214) * y -0.010557625006) * y
1066             +0.011630447319) * y -0.009279453341) * y
1067           +0.005353579108) * y -0.002141268741) * y
1068         +0.000535310849) * y +0.999936657524
1069    if z > 0.0:
1070    prob = ((x+1.0)*0.5)
1071    else:
1072    prob = ((1.0-x)*0.5)
1073    return prob
1074
1075   
Note: See TracBrowser for help on using the repository browser.