Revision 10346:c99dada8a093, 39.4 KB checked in by janezd <janez.demsar@…>, 2 years ago (diff)

Polished documentation for logistic regression

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