source: orange/Orange/classification/logreg.py @ 10246:11b418321f79

Revision 10246:11b418321f79, 39.8 KB checked in by janezd <janez.demsar@…>, 2 years ago (diff)

Unified argument names in and 2.5 and 3.0; numerous other changes in documentation

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 of all major features in logistic regression
11    classifier.
12
13    :param classifier: logistic regression classifier.
14    """
15
16    # print out class values
17    out = ['']
18    out.append("class attribute = " + classifier.domain.class_var.name)
19    out.append("class values = " + str(classifier.domain.class_var.values))
20    out.append('')
21   
22    # get the longest attribute name
23    longest=0
24    for at in classifier.continuized_domain.features:
25        if len(at.name)>longest:
26            longest=len(at.name)
27
28    # print out the head
29    formatstr = "%"+str(longest)+"s %10s %10s %10s %10s %10s"
30    out.append(formatstr % ("Feature", "beta", "st. error", "wald Z", "P", "OR=exp(beta)"))
31    out.append('')
32    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f"   
33    out.append(formatstr % ("Intercept", classifier.beta[0], classifier.beta_se[0], classifier.wald_Z[0], classifier.P[0]))
34    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f %10.2f"   
35    for i in range(len(classifier.continuized_domain.features)):
36        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])))
37
38    return '\n'.join(out)
39       
40
41def has_discrete_values(domain):
42    """
43    Return 1 if the given domain contains any discrete features, else 0.
44
45    :param domain: domain.
46    :type domain: :class:`Orange.data.Domain`
47    """
48    return any(at.var_type == Orange.feature.Type.Discrete
49               for at in domain.features)
50
51
52class LogRegLearner(Orange.classification.Learner):
53    """ Logistic regression learner.
54
55    If data instances are provided to
56    the constructor, the learning algorithm is called and the resulting
57    classifier is returned instead of the learner.
58
59    :param data: data table with either discrete or continuous features
60    :type data: Orange.data.Table
61    :param weight_id: the ID of the weight meta attribute
62    :type weight_id: int
63    :param remove_singular: set to 1 if you want automatic removal of
64        disturbing features, such as constants and singularities
65    :type remove_singular: bool
66    :param fitter: the fitting algorithm (by default the Newton-Raphson
67        fitting algorithm is used)
68    :param stepwise_lr: set to 1 if you wish to use stepwise logistic
69        regression
70    :type stepwise_lr: bool
71    :param add_crit: parameter for stepwise feature selection
72    :type add_crit: float
73    :param delete_crit: parameter for stepwise feature selection
74    :type delete_crit: float
75    :param num_features: parameter for stepwise feature selection
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        """Learn from the given table of data instances.
99
100        :param data: Data instances to learn from.
101        :type data: :class:`~Orange.data.Table`
102        :param weight: Id of meta attribute with weights of instances
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):
114            add_crit = getattr(self, "add_crit", 0.2)
115            delete_crit = getattr(self, "delete_crit", 0.3)
116            num_features = getattr(self, "num_features", -1)
117            attributes = StepWiseFSS(data, add_crit= add_crit,
118                delete_crit=delete_crit, imputer = imputer, num_features= num_features)
119            tmp_domain = Orange.data.Domain(attributes,
120                data.domain.class_var)
121            tmp_domain.addmetas(data.domain.getmetas())
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)
144            new_domain.addmetas(data.domain.getmetas())
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",
152                                    "addCrit": "add_crit",
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])
221                    new_domain.addmetas(data.domain.getmetas())
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])
235                    new_domain.addmetas(data.domain.getmetas())
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()
258            data.addMetaAttribute(weight, 1.0)
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:
278            beta_add = m.beta[m.continuized_domain.features[-1]]
279            betas_ap.append(beta_add)
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])
341                    newDomain.addmetas(newData.domain.getmetas())
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])
358                    newDomain.addmetas(orig_data.domain.getmetas())
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()
381            data.addMetaAttribute(weight, 1.0)
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:
396            beta_add = m.beta[m.continuized_domain.features[-1]]
397            betas_ap.append(beta_add)
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  Algorithm described in Hosmer and Lemeshow,
688  Applied Logistic Regression, 2000.
689
690  Perform stepwise logistic regression and return a list of the
691  most "informative" features. Each step of the algorithm is composed
692  of two parts. The first is backward elimination, where each already
693  chosen feature is tested for a significant contribution to the overall
694  model. If the worst among all tested features has higher significance
695  than is specified in :obj:`delete_crit`, the feature is removed from
696  the model. The second step is forward selection, which is similar to
697  backward elimination. It loops through all the features that are not
698  in the model and tests whether they contribute to the common model
699  with significance lower that :obj:`add_crit`. The algorithm stops when
700  no feature in the model is to be removed and no feature not in the
701  model is to be added. By setting :obj:`num_features` larger than -1,
702  the algorithm will stop its execution when the number of features in model
703  exceeds that number.
704
705  Significances are assesed via the likelihood ration chi-square
706  test. Normal F test is not appropriate, because errors are assumed to
707  follow a binomial distribution.
708
709  If :obj:`table` is specified, stepwise logistic regression implemented
710  in :obj:`StepWiseFSS` is performed and a list of chosen features
711  is returned. If :obj:`table` is not specified, an instance of
712  :obj:`StepWiseFSS` with all parameters set is returned and can be called
713  with data later.
714
715  :param table: data set.
716  :type table: Orange.data.Table
717
718  :param add_crit: "Alpha" level to judge if variable has enough importance to
719       be added in the new set. (e.g. if add_crit is 0.2,
720       then features is added if its P is lower than 0.2).
721  :type add_crit: float
722
723  :param delete_crit: Similar to add_crit, just that it is used at backward
724      elimination. It should be higher than add_crit!
725  :type delete_crit: float
726
727  :param num_features: maximum number of selected features,
728      use -1 for infinity.
729  :type num_features: int
730  :rtype: :obj:`StepWiseFSS` or list of features
731
732  """
733
734  def __new__(cls, data=None, **argkw):
735      self = Orange.classification.Learner.__new__(cls, **argkw)
736      if data:
737          self.__init__(**argkw)
738          return self.__call__(data)
739      else:
740          return self
741
742  @deprecated_keywords({"addCrit": "add_crit", "deleteCrit": "delete_crit",
743                        "numFeatures": "num_features"})
744  def __init__(self, add_crit=0.2, delete_crit=0.3, num_features = -1, **kwds):
745    self.__dict__.update(kwds)
746    self.add_crit = add_crit
747    self.delete_crit = delete_crit
748    self.num_features = num_features
749
750  def __call__(self, examples):
751    if getattr(self, "imputer", 0):
752        examples = self.imputer(examples)(examples)
753    if getattr(self, "removeMissing", 0):
754        examples = Orange.core.Preprocessor_dropMissing(examples)
755    continuizer = Orange.preprocess.DomainContinuizer(zeroBased=1,
756        continuousTreatment=Orange.preprocess.DomainContinuizer.Leave,
757                                           multinomialTreatment = Orange.preprocess.DomainContinuizer.FrequentIsBase,
758                                           classTreatment = Orange.preprocess.DomainContinuizer.Ignore)
759    attr = []
760    remain_attr = examples.domain.features[:]
761
762    # get LL for Majority Learner
763    tempDomain = Orange.data.Domain(attr,examples.domain.class_var)
764    #tempData  = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
765    tempData  = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
766
767    ll_Old = get_likelihood(LogRegFitter_Cholesky(), tempData)
768    ll_Best = -1000000
769    length_Old = float(len(tempData))
770
771    stop = 0
772    while not stop:
773        # LOOP until all variables are added or no further deletion nor addition of attribute is possible
774        worstAt = None
775        # if there are more than 1 attribute then perform backward elimination
776        if len(attr) >= 2:
777            minG = 1000
778            worstAt = attr[0]
779            ll_Best = ll_Old
780            length_Best = length_Old
781            for at in attr:
782                # check all attribute whether its presence enough increases LL?
783
784                tempAttr = filter(lambda x: x!=at, attr)
785                tempDomain = Orange.data.Domain(tempAttr,examples.domain.class_var)
786                tempDomain.addmetas(examples.domain.getmetas())
787                # domain, calculate P for LL improvement.
788                tempDomain  = continuizer(Orange.core.Preprocessor_dropMissing(examples.select(tempDomain)))
789                tempData = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
790
791                ll_Delete = get_likelihood(LogRegFitter_Cholesky(), tempData)
792                length_Delete = float(len(tempData))
793                length_Avg = (length_Delete + length_Old)/2.0
794
795                G=-2*length_Avg*(ll_Delete/length_Delete-ll_Old/length_Old)
796
797                # set new worst attribute
798                if G<minG:
799                    worstAt = at
800                    minG=G
801                    ll_Best = ll_Delete
802                    length_Best = length_Delete
803            # deletion of attribute
804
805            if worstAt.var_type==Orange.feature.Type.Continuous:
806                P=lchisqprob(minG,1);
807            else:
808                P=lchisqprob(minG,len(worstAt.values)-1);
809            if P>=self.delete_crit:
810                attr.remove(worstAt)
811                remain_attr.append(worstAt)
812                nodeletion=0
813                ll_Old = ll_Best
814                length_Old = length_Best
815            else:
816                nodeletion=1
817        else:
818            nodeletion = 1
819            # END OF DELETION PART
820
821        # if enough attributes has been chosen, stop the procedure
822        if self.num_features>-1 and len(attr)>=self.num_features:
823            remain_attr=[]
824
825        # for each attribute in the remaining
826        maxG=-1
827        ll_Best = ll_Old
828        length_Best = length_Old
829        bestAt = None
830        for at in remain_attr:
831            tempAttr = attr + [at]
832            tempDomain = Orange.data.Domain(tempAttr,examples.domain.class_var)
833            tempDomain.addmetas(examples.domain.getmetas())
834            # domain, calculate P for LL improvement.
835            tempDomain  = continuizer(Orange.core.Preprocessor_dropMissing(examples.select(tempDomain)))
836            tempData = Orange.core.Preprocessor_dropMissing(examples.select(tempDomain))
837            ll_New = get_likelihood(LogRegFitter_Cholesky(), tempData)
838
839            length_New = float(len(tempData)) # get number of examples in tempData to normalize likelihood
840
841            # P=PR(CHI^2>G), G=-2(L(0)-L(1))=2(E(0)-E(1))
842            length_avg = (length_New + length_Old)/2
843            G=-2*length_avg*(ll_Old/length_Old-ll_New/length_New);
844            if G>maxG:
845                bestAt = at
846                maxG=G
847                ll_Best = ll_New
848                length_Best = length_New
849        if not bestAt:
850            stop = 1
851            continue
852
853        if bestAt.var_type==Orange.feature.Type.Continuous:
854            P=lchisqprob(maxG,1);
855        else:
856            P=lchisqprob(maxG,len(bestAt.values)-1);
857        # Add attribute with smallest P to attributes(attr)
858        if P<=self.add_crit:
859            attr.append(bestAt)
860            remain_attr.remove(bestAt)
861            ll_Old = ll_Best
862            length_Old = length_Best
863
864        if (P>self.add_crit and nodeletion) or (bestAt == worstAt):
865            stop = 1
866
867    return attr
868
869StepWiseFSS = deprecated_members({"addCrit": "add_crit",
870                                   "deleteCrit": "delete_crit",
871                                   "numFeatures": "num_features"})(StepWiseFSS)
872
873
874class StepWiseFSSFilter(object):
875    def __new__(cls, data=None, **argkw):
876        self = object.__new__(cls)
877        if data:
878            self.__init__(**argkw)
879            return self.__call__(data)
880        else:
881            return self
882
883    @deprecated_keywords({"addCrit": "add_crit", "deleteCrit": "delete_crit",
884                          "numFeatures": "num_features"})
885    def __init__(self, add_crit=0.2, delete_crit=0.3, num_features = -1):
886        self.add_crit = add_crit
887        self.delete_crit = delete_crit
888        self.num_features = num_features
889
890    @deprecated_keywords({"examples": "data"})
891    def __call__(self, data):
892        attr = StepWiseFSS(data, add_crit=self.add_crit,
893            delete_crit= self.delete_crit, num_features= self.num_features)
894        return data.select(Orange.data.Domain(attr, data.domain.class_var))
895
896StepWiseFSSFilter = deprecated_members({"addCrit": "add_crit",
897                                        "deleteCrit": "delete_crit",
898                                        "numFeatures": "num_features"})\
899    (StepWiseFSSFilter)
900
901
902####################################
903##  PROBABILITY CALCULATIONS
904
905def lchisqprob(chisq,df):
906    """
907    Return the (1-tailed) probability value associated with the provided
908    chi-square value and df.  Adapted from chisq.c in Gary Perlman's |Stat.
909    """
910    BIG = 20.0
911    def ex(x):
912        BIG = 20.0
913        if x < -BIG:
914            return 0.0
915        else:
916            return math.exp(x)
917    if chisq <=0 or df < 1:
918        return 1.0
919    a = 0.5 * chisq
920    if df%2 == 0:
921        even = 1
922    else:
923        even = 0
924    if df > 1:
925        y = ex(-a)
926    if even:
927        s = y
928    else:
929        s = 2.0 * zprob(-math.sqrt(chisq))
930    if (df > 2):
931        chisq = 0.5 * (df - 1.0)
932        if even:
933            z = 1.0
934        else:
935            z = 0.5
936        if a > BIG:
937            if even:
938                e = 0.0
939            else:
940                e = math.log(math.sqrt(math.pi))
941            c = math.log(a)
942            while (z <= chisq):
943                e = math.log(z) + e
944                s = s + ex(c*z-a-e)
945                z = z + 1.0
946            return s
947        else:
948            if even:
949                e = 1.0
950            else:
951                e = 1.0 / math.sqrt(math.pi) / math.sqrt(a)
952            c = 0.0
953            while (z <= chisq):
954                e = e * (a/float(z))
955                c = c + e
956                z = z + 1.0
957            return (c*y+s)
958    else:
959        return s
960
961
962def zprob(z):
963    """
964    Returns the area under the normal curve 'to the left of' the given z value.
965    Thus::
966
967    for z<0, zprob(z) = 1-tail probability
968    for z>0, 1.0-zprob(z) = 1-tail probability
969    for any z, 2.0*(1.0-zprob(abs(z))) = 2-tail probability
970
971    Adapted from z.c in Gary Perlman's |Stat.
972    """
973    Z_MAX = 6.0    # maximum meaningful z-value
974    if z == 0.0:
975    x = 0.0
976    else:
977    y = 0.5 * math.fabs(z)
978    if y >= (Z_MAX*0.5):
979        x = 1.0
980    elif (y < 1.0):
981        w = y*y
982        x = ((((((((0.000124818987 * w
983            -0.001075204047) * w +0.005198775019) * w
984              -0.019198292004) * w +0.059054035642) * w
985            -0.151968751364) * w +0.319152932694) * w
986          -0.531923007300) * w +0.797884560593) * y * 2.0
987    else:
988        y = y - 2.0
989        x = (((((((((((((-0.000045255659 * y
990                 +0.000152529290) * y -0.000019538132) * y
991               -0.000676904986) * y +0.001390604284) * y
992             -0.000794620820) * y -0.002034254874) * y
993               +0.006549791214) * y -0.010557625006) * y
994             +0.011630447319) * y -0.009279453341) * y
995           +0.005353579108) * y -0.002141268741) * y
996         +0.000535310849) * y +0.999936657524
997    if z > 0.0:
998    prob = ((x+1.0)*0.5)
999    else:
1000    prob = ((1.0-x)*0.5)
1001    return prob
1002
1003   
Note: See TracBrowser for help on using the repository browser.