# source:orange/Orange/classification/logreg.py@11063:68bc42e33a2a

Revision 11063:68bc42e33a2a, 42.2 KB checked in by astaric, 16 months ago (diff)

Python 2.6 compatibility.

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