source: orange/Orange/classification/logreg.py @ 10883:79fbadeec704

Revision 10883:79fbadeec704, 41.5 KB checked in by markotoplak, 23 months ago (diff)

Logistic regression with LibLinear: added an explicit bias parameter, some documentation changes.

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