source: orange/Orange/classification/logreg.py @ 9878:80400e49d3a4

Revision 9878:80400e49d3a4, 40.1 KB checked in by Matija Polajnar <matija.polajnar@…>, 2 years ago (diff)

Remove arguments from object.new calls. It gives a deprecation warning and has no effects: the arguments get passed anyway.

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