source: orange/Orange/classification/logreg.py @ 11060:340b8bf1cbb4

Revision 11060:340b8bf1cbb4, 42.0 KB checked in by Yuval Greenfield <yuval@…>, 16 months ago (diff)

Avoid OverflowError when the beta is really high

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