Changeset 7193:1206ec12c326 in orange


Ignore:
Timestamp:
02/02/11 16:12:38 (3 years ago)
Author:
jzbontar <jure.zbontar@…>
Branch:
default
Convert:
754aa8ea13940bfa2756163980288cf1ed5e227e
Message:

logreg renaming init

Location:
orange
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • orange/Orange/classify/logreg.py

    r6848 r7193  
     1""" 
     2.. index: logreg 
     3 
     4=================== 
     5Logistic Regression 
     6=================== 
     7 
     8Module orngLR is a set of wrappers around classes LogisticLearner and 
     9LogisticClassifier, that are implemented in core Orange. This module 
     10expanses use of logistic regression to discrete attributes, it helps 
     11handling various anomalies in attributes, such as constant variables 
     12and singularities, that makes fitting logistic regression almost 
     13impossible. It also implements a function for constructing a stepwise 
     14logistic regression, which is a good technique to prevent overfitting, 
     15and is a good feature subset selection technique as well. 
     16 
     17""" 
     18 
    119from orange import LogRegLearner, LogRegClassifier, LogRegFitter, LogRegFitter_Cholesky 
     20 
     21import orange 
     22import orngCI 
     23import math, os 
     24from numpy import * 
     25from numpy.linalg import * 
     26 
     27 
     28####################### 
     29## Print out methods ## 
     30####################### 
     31 
     32def printOUT(classifier): 
     33    # print out class values 
     34    print 
     35    print "class attribute = " + classifier.domain.classVar.name 
     36    print "class values = " + str(classifier.domain.classVar.values) 
     37    print 
     38     
     39    # get the longest attribute name 
     40    longest=0 
     41    for at in classifier.continuizedDomain.attributes: 
     42        if len(at.name)>longest: 
     43            longest=len(at.name); 
     44 
     45    # print out the head 
     46    formatstr = "%"+str(longest)+"s %10s %10s %10s %10s %10s" 
     47    print formatstr % ("Attribute", "beta", "st. error", "wald Z", "P", "OR=exp(beta)") 
     48    print 
     49    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f"     
     50    print formatstr % ("Intercept", classifier.beta[0], classifier.beta_se[0], classifier.wald_Z[0], classifier.P[0]) 
     51    formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f %10.2f"     
     52    for i in range(len(classifier.continuizedDomain.attributes)): 
     53        print formatstr % (classifier.continuizedDomain.attributes[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])) 
     54         
     55 
     56def hasDiscreteValues(domain): 
     57    for at in domain.attributes: 
     58        if at.varType == orange.VarTypes.Discrete: 
     59            return 1 
     60    return 0 
     61 
     62def LogRegLearner(examples = None, weightID=0, **kwds): 
     63    lr = LogRegLearnerClass(**kwds) 
     64    if examples: 
     65        return lr(examples, weightID) 
     66    else: 
     67        return lr 
     68 
     69class LogRegLearnerClass(orange.Learner): 
     70    def __init__(self, removeSingular=0, fitter = None, **kwds): 
     71        self.__dict__.update(kwds) 
     72        self.removeSingular = removeSingular 
     73        self.fitter = None 
     74 
     75    def __call__(self, examples, weight=0): 
     76        imputer = getattr(self, "imputer", None) or None 
     77        if getattr(self, "removeMissing", 0): 
     78            examples = orange.Preprocessor_dropMissing(examples) 
     79##        if hasDiscreteValues(examples.domain): 
     80##            examples = createNoDiscTable(examples) 
     81        if not len(examples): 
     82            return None 
     83        if getattr(self, "stepwiseLR", 0): 
     84            addCrit = getattr(self, "addCrit", 0.2) 
     85            removeCrit = getattr(self, "removeCrit", 0.3) 
     86            numAttr = getattr(self, "numAttr", -1) 
     87            attributes = StepWiseFSS(examples, addCrit = addCrit, deleteCrit = removeCrit, imputer = imputer, numAttr = numAttr) 
     88            tmpDomain = orange.Domain(attributes, examples.domain.classVar) 
     89            tmpDomain.addmetas(examples.domain.getmetas()) 
     90            examples = examples.select(tmpDomain) 
     91        learner = orange.LogRegLearner() 
     92        learner.imputerConstructor = imputer 
     93        if imputer: 
     94            examples = self.imputer(examples)(examples) 
     95        examples = orange.Preprocessor_dropMissing(examples) 
     96        if self.fitter: 
     97            learner.fitter = self.fitter 
     98        if self.removeSingular: 
     99            lr = learner.fitModel(examples, weight) 
     100        else: 
     101            lr = learner(examples, weight) 
     102        while isinstance(lr, orange.Variable): 
     103            if isinstance(lr.getValueFrom, orange.ClassifierFromVar) and isinstance(lr.getValueFrom.transformer, orange.Discrete2Continuous): 
     104                lr = lr.getValueFrom.variable 
     105            attributes = examples.domain.attributes[:] 
     106            if lr in attributes: 
     107                attributes.remove(lr) 
     108            else: 
     109                attributes.remove(lr.getValueFrom.variable) 
     110            newDomain = orange.Domain(attributes, examples.domain.classVar) 
     111            newDomain.addmetas(examples.domain.getmetas()) 
     112            examples = examples.select(newDomain) 
     113            lr = learner.fitModel(examples, weight) 
     114        return lr 
     115 
     116 
     117 
     118def Univariate_LogRegLearner(examples=None, **kwds): 
     119    learner = apply(Univariate_LogRegLearner_Class, (), kwds) 
     120    if examples: 
     121        return learner(examples) 
     122    else: 
     123        return learner 
     124 
     125class Univariate_LogRegLearner_Class(orange.Learner): 
     126    def __init__(self, **kwds): 
     127        self.__dict__.update(kwds) 
     128 
     129    def __call__(self, examples): 
     130        examples = createFullNoDiscTable(examples) 
     131        classifiers = map(lambda x: LogRegLearner(orange.Preprocessor_dropMissing(examples.select(orange.Domain(x, examples.domain.classVar)))), examples.domain.attributes) 
     132        maj_classifier = LogRegLearner(orange.Preprocessor_dropMissing(examples.select(orange.Domain(examples.domain.classVar)))) 
     133        beta = [maj_classifier.beta[0]] + [x.beta[1] for x in classifiers] 
     134        beta_se = [maj_classifier.beta_se[0]] + [x.beta_se[1] for x in classifiers] 
     135        P = [maj_classifier.P[0]] + [x.P[1] for x in classifiers] 
     136        wald_Z = [maj_classifier.wald_Z[0]] + [x.wald_Z[1] for x in classifiers] 
     137        domain = examples.domain 
     138 
     139        return Univariate_LogRegClassifier(beta = beta, beta_se = beta_se, P = P, wald_Z = wald_Z, domain = domain) 
     140 
     141class Univariate_LogRegClassifier(orange.Classifier): 
     142    def __init__(self, **kwds): 
     143        self.__dict__.update(kwds) 
     144 
     145    def __call__(self, example, resultType = orange.GetValue): 
     146        # classification not implemented yet. For now its use is only to provide regression coefficients and its statistics 
     147        pass 
     148     
     149 
     150def LogRegLearner_getPriors(examples = None, weightID=0, **kwds): 
     151    lr = LogRegLearnerClass_getPriors(**kwds) 
     152    if examples: 
     153        return lr(examples, weightID) 
     154    else: 
     155        return lr 
     156 
     157class LogRegLearnerClass_getPriors(object): 
     158    def __init__(self, removeSingular=0, **kwds): 
     159        self.__dict__.update(kwds) 
     160        self.removeSingular = removeSingular 
     161    def __call__(self, examples, weight=0): 
     162        # next function changes data set to a extended with unknown values  
     163        def createLogRegExampleTable(data, weightID): 
     164            setsOfData = [] 
     165            for at in data.domain.attributes: 
     166                # za vsak atribut kreiraj nov newExampleTable newData 
     167                # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable 
     168                if at.varType == orange.VarTypes.Continuous: 
     169                    atDisc = orange.FloatVariable(at.name + "Disc") 
     170                    newDomain = orange.Domain(data.domain.attributes+[atDisc,data.domain.classVar]) 
     171                    newDomain.addmetas(data.domain.getmetas()) 
     172                    newData = orange.ExampleTable(newDomain,data) 
     173                    altData = orange.ExampleTable(newDomain,data) 
     174                    for i,d in enumerate(newData): 
     175                        d[atDisc] = 0 
     176                        d[weightID] = 1*data[i][weightID] 
     177                    for i,d in enumerate(altData): 
     178                        d[atDisc] = 1 
     179                        d[at] = 0 
     180                        d[weightID] = 0.000001*data[i][weightID] 
     181                elif at.varType == orange.VarTypes.Discrete: 
     182                # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X" 
     183                    atNew = orange.EnumVariable(at.name, values = at.values + [at.name+"X"]) 
     184                    newDomain = orange.Domain(filter(lambda x: x!=at, data.domain.attributes)+[atNew,data.domain.classVar]) 
     185                    newDomain.addmetas(data.domain.getmetas()) 
     186                    newData = orange.ExampleTable(newDomain,data) 
     187                    altData = orange.ExampleTable(newDomain,data) 
     188                    for i,d in enumerate(newData): 
     189                        d[atNew] = data[i][at] 
     190                        d[weightID] = 1*data[i][weightID] 
     191                    for i,d in enumerate(altData): 
     192                        d[atNew] = at.name+"X" 
     193                        d[weightID] = 0.000001*data[i][weightID] 
     194                newData.extend(altData) 
     195                setsOfData.append(newData) 
     196            return setsOfData 
     197                   
     198        learner = LogRegLearner(imputer = orange.ImputerConstructor_average(), removeSingular = self.removeSingular) 
     199        # get Original Model 
     200        orig_model = learner(examples,weight) 
     201        if orig_model.fit_status: 
     202            print "Warning: model did not converge" 
     203 
     204        # get extended Model (you should not change data) 
     205        if weight == 0: 
     206            weight = orange.newmetaid() 
     207            examples.addMetaAttribute(weight, 1.0) 
     208        extended_set_of_examples = createLogRegExampleTable(examples, weight) 
     209        extended_models = [learner(extended_examples, weight) \ 
     210                           for extended_examples in extended_set_of_examples] 
     211 
     212##        print examples[0] 
     213##        printOUT(orig_model) 
     214##        print orig_model.domain 
     215##        print orig_model.beta 
     216##        print orig_model.beta[orig_model.continuizedDomain.attributes[-1]] 
     217##        for i,m in enumerate(extended_models): 
     218##            print examples.domain.attributes[i] 
     219##            printOUT(m) 
     220             
     221         
     222        # izracunas odstopanja 
     223        # get sum of all betas 
     224        beta = 0 
     225        betas_ap = [] 
     226        for m in extended_models: 
     227            beta_add = m.beta[m.continuizedDomain.attributes[-1]] 
     228            betas_ap.append(beta_add) 
     229            beta = beta + beta_add 
     230         
     231        # substract it from intercept 
     232        #print "beta", beta 
     233        logistic_prior = orig_model.beta[0]+beta 
     234         
     235        # compare it to bayes prior 
     236        bayes = orange.BayesLearner(examples) 
     237        bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0]) 
     238 
     239        # normalize errors 
     240##        print "bayes", bayes_prior 
     241##        print "lr", orig_model.beta[0] 
     242##        print "lr2", logistic_prior 
     243##        print "dist", orange.Distribution(examples.domain.classVar,examples) 
     244##        print "prej", betas_ap 
     245 
     246        # error normalization - to avoid errors due to assumption of independence of unknown values 
     247        dif = bayes_prior - logistic_prior 
     248        positives = sum(filter(lambda x: x>=0, betas_ap)) 
     249        negatives = -sum(filter(lambda x: x<0, betas_ap)) 
     250        if not negatives == 0: 
     251            kPN = positives/negatives 
     252            diffNegatives = dif/(1+kPN) 
     253            diffPositives = kPN*diffNegatives 
     254            kNegatives = (negatives-diffNegatives)/negatives 
     255            kPositives = positives/(positives-diffPositives) 
     256    ##        print kNegatives 
     257    ##        print kPositives 
     258 
     259            for i,b in enumerate(betas_ap): 
     260                if b<0: betas_ap[i]*=kNegatives 
     261                else: betas_ap[i]*=kPositives 
     262        #print "potem", betas_ap 
     263 
     264        # vrni originalni model in pripadajoce apriorne niclele 
     265        return (orig_model, betas_ap) 
     266        #return (bayes_prior,orig_model.beta[examples.domain.classVar],logistic_prior) 
     267 
     268class LogRegLearnerClass_getPriors_OneTable: 
     269    def __init__(self, removeSingular=0, **kwds): 
     270        self.__dict__.update(kwds) 
     271        self.removeSingular = removeSingular 
     272    def __call__(self, examples, weight=0): 
     273        # next function changes data set to a extended with unknown values  
     274        def createLogRegExampleTable(data, weightID): 
     275            finalData = orange.ExampleTable(data) 
     276            origData = orange.ExampleTable(data) 
     277            for at in data.domain.attributes: 
     278                # za vsak atribut kreiraj nov newExampleTable newData 
     279                # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable 
     280                if at.varType == orange.VarTypes.Continuous: 
     281                    atDisc = orange.FloatVariable(at.name + "Disc") 
     282                    newDomain = orange.Domain(origData.domain.attributes+[atDisc,data.domain.classVar]) 
     283                    newDomain.addmetas(newData.domain.getmetas()) 
     284                    finalData = orange.ExampleTable(newDomain,finalData) 
     285                    newData = orange.ExampleTable(newDomain,origData) 
     286                    origData = orange.ExampleTable(newDomain,origData) 
     287                    for d in origData: 
     288                        d[atDisc] = 0 
     289                    for d in finalData: 
     290                        d[atDisc] = 0 
     291                    for i,d in enumerate(newData): 
     292                        d[atDisc] = 1 
     293                        d[at] = 0 
     294                        d[weightID] = 100*data[i][weightID] 
     295                         
     296                elif at.varType == orange.VarTypes.Discrete: 
     297                # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X" 
     298                    atNew = orange.EnumVariable(at.name, values = at.values + [at.name+"X"]) 
     299                    newDomain = orange.Domain(filter(lambda x: x!=at, origData.domain.attributes)+[atNew,origData.domain.classVar]) 
     300                    newDomain.addmetas(origData.domain.getmetas()) 
     301                    temp_finalData = orange.ExampleTable(finalData) 
     302                    finalData = orange.ExampleTable(newDomain,finalData) 
     303                    newData = orange.ExampleTable(newDomain,origData) 
     304                    temp_origData = orange.ExampleTable(origData) 
     305                    origData = orange.ExampleTable(newDomain,origData) 
     306                    for i,d in enumerate(origData): 
     307                        d[atNew] = temp_origData[i][at] 
     308                    for i,d in enumerate(finalData): 
     309                        d[atNew] = temp_finalData[i][at]                         
     310                    for i,d in enumerate(newData): 
     311                        d[atNew] = at.name+"X" 
     312                        d[weightID] = 10*data[i][weightID] 
     313                finalData.extend(newData) 
     314            return finalData 
     315                   
     316        learner = LogRegLearner(imputer = orange.ImputerConstructor_average(), removeSingular = self.removeSingular) 
     317        # get Original Model 
     318        orig_model = learner(examples,weight) 
     319 
     320        # get extended Model (you should not change data) 
     321        if weight == 0: 
     322            weight = orange.newmetaid() 
     323            examples.addMetaAttribute(weight, 1.0) 
     324        extended_examples = createLogRegExampleTable(examples, weight) 
     325        extended_model = learner(extended_examples, weight) 
     326 
     327##        print examples[0] 
     328##        printOUT(orig_model) 
     329##        print orig_model.domain 
     330##        print orig_model.beta 
     331 
     332##        printOUT(extended_model)         
     333        # izracunas odstopanja 
     334        # get sum of all betas 
     335        beta = 0 
     336        betas_ap = [] 
     337        for m in extended_models: 
     338            beta_add = m.beta[m.continuizedDomain.attributes[-1]] 
     339            betas_ap.append(beta_add) 
     340            beta = beta + beta_add 
     341         
     342        # substract it from intercept 
     343        #print "beta", beta 
     344        logistic_prior = orig_model.beta[0]+beta 
     345         
     346        # compare it to bayes prior 
     347        bayes = orange.BayesLearner(examples) 
     348        bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0]) 
     349 
     350        # normalize errors 
     351        #print "bayes", bayes_prior 
     352        #print "lr", orig_model.beta[0] 
     353        #print "lr2", logistic_prior 
     354        #print "dist", orange.Distribution(examples.domain.classVar,examples) 
     355        k = (bayes_prior-orig_model.beta[0])/(logistic_prior-orig_model.beta[0]) 
     356        #print "prej", betas_ap 
     357        betas_ap = [k*x for x in betas_ap]                 
     358        #print "potem", betas_ap 
     359 
     360        # vrni originalni model in pripadajoce apriorne niclele 
     361        return (orig_model, betas_ap) 
     362        #return (bayes_prior,orig_model.beta[data.domain.classVar],logistic_prior) 
     363 
     364 
     365###################################### 
     366#### Fitters for logistic regression (logreg) learner #### 
     367###################################### 
     368 
     369def Pr(x, betas): 
     370    k = math.exp(dot(x, betas)) 
     371    return k / (1+k) 
     372 
     373def lh(x,y,betas): 
     374    llh = 0.0 
     375    for i,x_i in enumerate(x): 
     376        pr = Pr(x_i,betas) 
     377        llh += y[i]*log(max(pr,1e-6)) + (1-y[i])*log(max(1-pr,1e-6)) 
     378    return llh 
     379 
     380 
     381def diag(vector): 
     382    mat = identity(len(vector), Float) 
     383    for i,v in enumerate(vector): 
     384        mat[i][i] = v 
     385    return mat 
     386     
     387class simpleFitter(orange.LogRegFitter): 
     388    def __init__(self, penalty=0, se_penalty = False): 
     389        self.penalty = penalty 
     390        self.se_penalty = se_penalty 
     391    def __call__(self, data, weight=0): 
     392        ml = data.native(0) 
     393        for i in range(len(data.domain.attributes)): 
     394          a = data.domain.attributes[i] 
     395          if a.varType == orange.VarTypes.Discrete: 
     396            for m in ml: 
     397              m[i] = a.values.index(m[i]) 
     398        for m in ml: 
     399          m[-1] = data.domain.classVar.values.index(m[-1]) 
     400        Xtmp = array(ml) 
     401        y = Xtmp[:,-1]   # true probabilities (1's or 0's) 
     402        one = reshape(array([1]*len(data)), (len(data),1)) # intercept column 
     403        X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data 
     404 
     405        betas = array([0.0] * (len(data.domain.attributes)+1)) 
     406        oldBetas = array([1.0] * (len(data.domain.attributes)+1)) 
     407        N = len(data) 
     408 
     409        pen_matrix = array([self.penalty] * (len(data.domain.attributes)+1)) 
     410        if self.se_penalty: 
     411            p = array([Pr(X[i], betas) for i in range(len(data))]) 
     412            W = identity(len(data), Float) 
     413            pp = p * (1.0-p) 
     414            for i in range(N): 
     415                W[i,i] = pp[i] 
     416            se = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))))) 
     417            for i,p in enumerate(pen_matrix): 
     418                pen_matrix[i] *= se[i] 
     419        # predict the probability for an instance, x and betas are vectors 
     420        # start the computation 
     421        likelihood = 0. 
     422        likelihood_new = 1. 
     423        while abs(likelihood - likelihood_new)>1e-5: 
     424            likelihood = likelihood_new 
     425            oldBetas = betas 
     426            p = array([Pr(X[i], betas) for i in range(len(data))]) 
     427 
     428            W = identity(len(data), Float) 
     429            pp = p * (1.0-p) 
     430            for i in range(N): 
     431                W[i,i] = pp[i] 
     432 
     433            WI = inverse(W) 
     434            z = matrixmultiply(X, betas) + matrixmultiply(WI, y - p) 
     435 
     436            tmpA = inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))+diag(pen_matrix)) 
     437            tmpB = matrixmultiply(transpose(X), y-p) 
     438            betas = oldBetas + matrixmultiply(tmpA,tmpB) 
     439#            betaTemp = matrixmultiply(matrixmultiply(matrixmultiply(matrixmultiply(tmpA,transpose(X)),W),X),oldBetas) 
     440#            print betaTemp 
     441#            tmpB = matrixmultiply(transpose(X), matrixmultiply(W, z)) 
     442#            betas = matrixmultiply(tmpA, tmpB) 
     443            likelihood_new = lh(X,y,betas)-self.penalty*sum([b*b for b in betas]) 
     444            print likelihood_new 
     445 
     446             
     447             
     448##        XX = sqrt(diagonal(inverse(matrixmultiply(transpose(X),X)))) 
     449##        yhat = array([Pr(X[i], betas) for i in range(len(data))]) 
     450##        ss = sum((y - yhat) ** 2) / (N - len(data.domain.attributes) - 1) 
     451##        sigma = math.sqrt(ss) 
     452        p = array([Pr(X[i], betas) for i in range(len(data))]) 
     453        W = identity(len(data), Float) 
     454        pp = p * (1.0-p) 
     455        for i in range(N): 
     456            W[i,i] = pp[i] 
     457        diXWX = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))))) 
     458        xTemp = matrixmultiply(matrixmultiply(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))),transpose(X)),y) 
     459        beta = [] 
     460        beta_se = [] 
     461        print "likelihood ridge", likelihood 
     462        for i in range(len(betas)): 
     463            beta.append(betas[i]) 
     464            beta_se.append(diXWX[i]) 
     465        return (self.OK, beta, beta_se, 0) 
     466 
     467def Pr_bx(bx): 
     468    if bx > 35: 
     469        return 1 
     470    if bx < -35: 
     471        return 0 
     472    return exp(bx)/(1+exp(bx)) 
     473 
     474class bayesianFitter(orange.LogRegFitter): 
     475    def __init__(self, penalty=0, anch_examples=[], tau = 0): 
     476        self.penalty = penalty 
     477        self.anch_examples = anch_examples 
     478        self.tau = tau 
     479 
     480    def createArrayData(self,data): 
     481        if not len(data): 
     482            return (array([]),array([])) 
     483        # convert data to numeric 
     484        ml = data.native(0) 
     485        for i,a in enumerate(data.domain.attributes): 
     486          if a.varType == orange.VarTypes.Discrete: 
     487            for m in ml: 
     488              m[i] = a.values.index(m[i]) 
     489        for m in ml: 
     490          m[-1] = data.domain.classVar.values.index(m[-1]) 
     491        Xtmp = array(ml) 
     492        y = Xtmp[:,-1]   # true probabilities (1's or 0's) 
     493        one = reshape(array([1]*len(data)), (len(data),1)) # intercept column 
     494        X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data 
     495        return (X,y) 
     496     
     497    def __call__(self, data, weight=0): 
     498        (X,y)=self.createArrayData(data) 
     499 
     500        exTable = orange.ExampleTable(data.domain) 
     501        for id,ex in self.anch_examples: 
     502            exTable.extend(orange.ExampleTable(ex,data.domain)) 
     503        (X_anch,y_anch)=self.createArrayData(exTable) 
     504 
     505        betas = array([0.0] * (len(data.domain.attributes)+1)) 
     506 
     507        likelihood,betas = self.estimateBeta(X,y,betas,[0]*(len(betas)),X_anch,y_anch) 
     508 
     509        # get attribute groups atGroup = [(startIndex, number of values), ...) 
     510        ats = data.domain.attributes 
     511        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:] 
     512        atGroup=[[0,0]] 
     513        for v_i,v in enumerate(atVec): 
     514            if v[1]==0: atGroup[-1][1]+=1 
     515            else:       atGroup.append([v_i,1]) 
     516         
     517        # compute zero values for attributes 
     518        sumB = 0. 
     519        for ag in atGroup: 
     520            X_temp = concatenate((X[:,:ag[0]+1],X[:,ag[0]+1+ag[1]:]),1) 
     521            if X_anch: 
     522                X_anch_temp = concatenate((X_anch[:,:ag[0]+1],X_anch[:,ag[0]+1+ag[1]:]),1) 
     523            else: X_anch_temp = X_anch 
     524##            print "1", concatenate((betas[:i+1],betas[i+2:])) 
     525##            print "2", betas 
     526            likelihood_temp,betas_temp=self.estimateBeta(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) 
     527            print "finBetas", betas, betas_temp 
     528            print "betas", betas[0], betas_temp[0] 
     529            sumB += betas[0]-betas_temp[0] 
     530        apriori = orange.Distribution(data.domain.classVar, data) 
     531        aprioriProb = apriori[0]/apriori.abs 
     532         
     533        print "koncni rezultat", sumB, math.log((1-aprioriProb)/aprioriProb), betas[0] 
     534             
     535        beta = [] 
     536        beta_se = [] 
     537        print "likelihood2", likelihood 
     538        for i in range(len(betas)): 
     539            beta.append(betas[i]) 
     540            beta_se.append(0.0) 
     541        return (self.OK, beta, beta_se, 0) 
     542 
     543      
     544         
     545    def estimateBeta(self,X,y,betas,const_betas,X_anch,y_anch): 
     546        N,N_anch = len(y),len(y_anch) 
     547        r,r_anch = array([dot(X[i], betas) for i in range(N)]),\ 
     548                   array([dot(X_anch[i], betas) for i in range(N_anch)]) 
     549        p    = array([Pr_bx(ri) for ri in r]) 
     550        X_sq = X*X 
     551 
     552        max_delta      = [1.]*len(const_betas) 
     553        likelihood     = -1.e+10 
     554        likelihood_new = -1.e+9 
     555        while abs(likelihood - likelihood_new)>0.01 and max(max_delta)>0.01: 
     556            likelihood = likelihood_new 
     557            print likelihood 
     558            betas_temp = [b for b in betas] 
     559            for j in range(len(betas)): 
     560                if const_betas[j]: continue 
     561                dl = matrixmultiply(X[:,j],transpose(y-p)) 
     562                for xi,x in enumerate(X_anch): 
     563                    dl += self.penalty*x[j]*(y_anch[xi] - Pr_bx(r_anch[xi]*self.penalty)) 
     564 
     565                ddl = matrixmultiply(X_sq[:,j],transpose(p*(1-p))) 
     566                for xi,x in enumerate(X_anch): 
     567                    ddl += self.penalty*x[j]*Pr_bx(r[xi]*self.penalty)*(1-Pr_bx(r[xi]*self.penalty)) 
     568 
     569                if j==0: 
     570                    dv = dl/max(ddl,1e-6) 
     571                elif betas[j] == 0: # special handling due to non-defined first and second derivatives 
     572                    dv = (dl-self.tau)/max(ddl,1e-6) 
     573                    if dv < 0: 
     574                        dv = (dl+self.tau)/max(ddl,1e-6) 
     575                        if dv > 0: 
     576                            dv = 0 
     577                else: 
     578                    dl -= sign(betas[j])*self.tau 
     579                    dv = dl/max(ddl,1e-6) 
     580                    if not sign(betas[j] + dv) == sign(betas[j]): 
     581                        dv = -betas[j] 
     582                dv = min(max(dv,-max_delta[j]),max_delta[j]) 
     583                r+= X[:,j]*dv 
     584                p = array([Pr_bx(ri) for ri in r]) 
     585                if N_anch: 
     586                    r_anch+=X_anch[:,j]*dv 
     587                betas[j] += dv 
     588                max_delta[j] = max(2*abs(dv),max_delta[j]/2) 
     589            likelihood_new = lh(X,y,betas) 
     590            for xi,x in enumerate(X_anch): 
     591                try: 
     592                    likelihood_new += y_anch[xi]*r_anch[xi]*self.penalty-log(1+exp(r_anch[xi]*self.penalty)) 
     593                except: 
     594                    likelihood_new += r_anch[xi]*self.penalty*(y_anch[xi]-1) 
     595            likelihood_new -= sum([abs(b) for b in betas[1:]])*self.tau 
     596            if likelihood_new < likelihood: 
     597                max_delta = [md/4 for md in max_delta] 
     598                likelihood_new = likelihood 
     599                likelihood = likelihood_new + 1. 
     600                betas = [b for b in betas_temp] 
     601        print "betas", betas 
     602        print "init_like", likelihood_new 
     603        print "pure_like", lh(X,y,betas) 
     604        return (likelihood,betas) 
     605     
     606############################################################ 
     607####  Feature subset selection for logistic regression  #### 
     608############################################################ 
     609 
     610 
     611def StepWiseFSS(examples = None, **kwds): 
     612    """ 
     613      Constructs and returns a new set of examples that includes a 
     614      class and attributes selected by stepwise logistic regression. This is an 
     615      implementation of algorithm described in [Hosmer and Lemeshow, Applied Logistic Regression, 2000] 
     616 
     617      examples: data set (ExampleTable)      
     618      addCrit: "Alpha" level to judge if variable has enough importance to be added in the new set. (e.g. if addCrit is 0.2, then attribute is added if its P is lower than 0.2) 
     619      deleteCrit: Similar to addCrit, just that it is used at backward elimination. It should be higher than addCrit! 
     620      numAttr: maximum number of selected attributes, use -1 for infinity 
     621    """ 
     622 
     623    fss = apply(StepWiseFSS_class, (), kwds) 
     624    if examples is not None: 
     625        return fss(examples) 
     626    else: 
     627        return fss 
     628 
     629def getLikelihood(fitter, examples): 
     630    res = fitter(examples) 
     631    if res[0] in [fitter.OK]: #, fitter.Infinity, fitter.Divergence]: 
     632       status, beta, beta_se, likelihood = res 
     633       if sum([abs(b) for b in beta])<sum([abs(b) for b in beta_se]): 
     634           return -100*len(examples) 
     635       return likelihood 
     636    else: 
     637       return -100*len(examples) 
     638         
     639 
     640 
     641class StepWiseFSS_class(orange.Learner): 
     642  def __init__(self, addCrit=0.2, deleteCrit=0.3, numAttr = -1, **kwds): 
     643    self.__dict__.update(kwds) 
     644    self.addCrit = addCrit 
     645    self.deleteCrit = deleteCrit 
     646    self.numAttr = numAttr 
     647  def __call__(self, examples): 
     648    if getattr(self, "imputer", 0): 
     649        examples = self.imputer(examples)(examples) 
     650    if getattr(self, "removeMissing", 0): 
     651        examples = orange.Preprocessor_dropMissing(examples) 
     652    continuizer = orange.DomainContinuizer(zeroBased=1,continuousTreatment=orange.DomainContinuizer.Leave, 
     653                                           multinomialTreatment = orange.DomainContinuizer.FrequentIsBase, 
     654                                           classTreatment = orange.DomainContinuizer.Ignore) 
     655    attr = [] 
     656    remain_attr = examples.domain.attributes[:] 
     657 
     658    # get LL for Majority Learner  
     659    tempDomain = orange.Domain(attr,examples.domain.classVar) 
     660    #tempData  = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
     661    tempData  = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
     662 
     663    ll_Old = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
     664    ll_Best = -1000000 
     665    length_Old = float(len(tempData)) 
     666 
     667    stop = 0 
     668    while not stop: 
     669        # LOOP until all variables are added or no further deletion nor addition of attribute is possible 
     670        worstAt = None 
     671        # if there are more than 1 attribute then perform backward elimination 
     672        if len(attr) >= 2: 
     673            minG = 1000 
     674            worstAt = attr[0] 
     675            ll_Best = ll_Old 
     676            length_Best = length_Old 
     677            for at in attr: 
     678                # check all attribute whether its presence enough increases LL? 
     679 
     680                tempAttr = filter(lambda x: x!=at, attr) 
     681                tempDomain = orange.Domain(tempAttr,examples.domain.classVar) 
     682                tempDomain.addmetas(examples.domain.getmetas()) 
     683                # domain, calculate P for LL improvement. 
     684                tempDomain  = continuizer(orange.Preprocessor_dropMissing(examples.select(tempDomain))) 
     685                tempData = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
     686 
     687                ll_Delete = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
     688                length_Delete = float(len(tempData)) 
     689                length_Avg = (length_Delete + length_Old)/2.0 
     690 
     691                G=-2*length_Avg*(ll_Delete/length_Delete-ll_Old/length_Old) 
     692 
     693                # set new worst attribute                 
     694                if G<minG: 
     695                    worstAt = at 
     696                    minG=G 
     697                    ll_Best = ll_Delete 
     698                    length_Best = length_Delete 
     699            # deletion of attribute 
     700             
     701            if worstAt.varType==orange.VarTypes.Continuous: 
     702                P=lchisqprob(minG,1); 
     703            else: 
     704                P=lchisqprob(minG,len(worstAt.values)-1); 
     705            if P>=self.deleteCrit: 
     706                attr.remove(worstAt) 
     707                remain_attr.append(worstAt) 
     708                nodeletion=0 
     709                ll_Old = ll_Best 
     710                length_Old = length_Best 
     711            else: 
     712                nodeletion=1 
     713        else: 
     714            nodeletion = 1 
     715            # END OF DELETION PART 
     716             
     717        # if enough attributes has been chosen, stop the procedure 
     718        if self.numAttr>-1 and len(attr)>=self.numAttr: 
     719            remain_attr=[] 
     720          
     721        # for each attribute in the remaining 
     722        maxG=-1 
     723        ll_Best = ll_Old 
     724        length_Best = length_Old 
     725        bestAt = None 
     726        for at in remain_attr: 
     727            tempAttr = attr + [at] 
     728            tempDomain = orange.Domain(tempAttr,examples.domain.classVar) 
     729            tempDomain.addmetas(examples.domain.getmetas()) 
     730            # domain, calculate P for LL improvement. 
     731            tempDomain  = continuizer(orange.Preprocessor_dropMissing(examples.select(tempDomain))) 
     732            tempData = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
     733            ll_New = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
     734 
     735            length_New = float(len(tempData)) # get number of examples in tempData to normalize likelihood 
     736 
     737            # P=PR(CHI^2>G), G=-2(L(0)-L(1))=2(E(0)-E(1)) 
     738            length_avg = (length_New + length_Old)/2 
     739            G=-2*length_avg*(ll_Old/length_Old-ll_New/length_New); 
     740            if G>maxG: 
     741                bestAt = at 
     742                maxG=G 
     743                ll_Best = ll_New 
     744                length_Best = length_New 
     745        if not bestAt: 
     746            stop = 1 
     747            continue 
     748         
     749        if bestAt.varType==orange.VarTypes.Continuous: 
     750            P=lchisqprob(maxG,1); 
     751        else: 
     752            P=lchisqprob(maxG,len(bestAt.values)-1); 
     753        # Add attribute with smallest P to attributes(attr) 
     754        if P<=self.addCrit: 
     755            attr.append(bestAt) 
     756            remain_attr.remove(bestAt) 
     757            ll_Old = ll_Best 
     758            length_Old = length_Best 
     759 
     760        if (P>self.addCrit and nodeletion) or (bestAt == worstAt): 
     761            stop = 1 
     762 
     763    return attr 
     764 
     765 
     766def StepWiseFSS_Filter(examples = None, **kwds): 
     767    """ 
     768        check function StepWiseFSS() 
     769    """ 
     770 
     771    filter = apply(StepWiseFSS_Filter_class, (), kwds) 
     772    if examples is not None: 
     773        return filter(examples) 
     774    else: 
     775        return filter 
     776 
     777 
     778class StepWiseFSS_Filter_class(object): 
     779    def __init__(self, addCrit=0.2, deleteCrit=0.3, numAttr = -1): 
     780        self.addCrit = addCrit 
     781        self.deleteCrit = deleteCrit 
     782        self.numAttr = numAttr 
     783    def __call__(self, examples): 
     784        attr = StepWiseFSS(examples, addCrit=self.addCrit, deleteCrit = self.deleteCrit, numAttr = self.numAttr) 
     785        return examples.select(orange.Domain(attr, examples.domain.classVar)) 
     786                 
     787 
     788#################################### 
     789####  PROBABILITY CALCULATIONS  #### 
     790#################################### 
     791 
     792def lchisqprob(chisq,df): 
     793    """ 
     794Returns the (1-tailed) probability value associated with the provided 
     795chi-square value and df.  Adapted from chisq.c in Gary Perlman's |Stat. 
     796 
     797Usage:   lchisqprob(chisq,df) 
     798""" 
     799    BIG = 20.0 
     800    def ex(x): 
     801        BIG = 20.0 
     802        if x < -BIG: 
     803            return 0.0 
     804        else: 
     805            return math.exp(x) 
     806    if chisq <=0 or df < 1: 
     807        return 1.0 
     808    a = 0.5 * chisq 
     809    if df%2 == 0: 
     810        even = 1 
     811    else: 
     812        even = 0 
     813    if df > 1: 
     814        y = ex(-a) 
     815    if even: 
     816        s = y 
     817    else: 
     818        s = 2.0 * zprob(-math.sqrt(chisq)) 
     819    if (df > 2): 
     820        chisq = 0.5 * (df - 1.0) 
     821        if even: 
     822            z = 1.0 
     823        else: 
     824            z = 0.5 
     825        if a > BIG: 
     826            if even: 
     827                e = 0.0 
     828            else: 
     829                e = math.log(math.sqrt(math.pi)) 
     830            c = math.log(a) 
     831            while (z <= chisq): 
     832                e = math.log(z) + e 
     833                s = s + ex(c*z-a-e) 
     834                z = z + 1.0 
     835            return s 
     836        else: 
     837            if even: 
     838                e = 1.0 
     839            else: 
     840                e = 1.0 / math.sqrt(math.pi) / math.sqrt(a) 
     841            c = 0.0 
     842            while (z <= chisq): 
     843                e = e * (a/float(z)) 
     844                c = c + e 
     845                z = z + 1.0 
     846            return (c*y+s) 
     847    else: 
     848        return s 
     849 
     850 
     851def zprob(z): 
     852    """ 
     853Returns the area under the normal curve 'to the left of' the given z value. 
     854Thus,  
     855    for z<0, zprob(z) = 1-tail probability 
     856    for z>0, 1.0-zprob(z) = 1-tail probability 
     857    for any z, 2.0*(1.0-zprob(abs(z))) = 2-tail probability 
     858Adapted from z.c in Gary Perlman's |Stat. 
     859 
     860Usage:   lzprob(z) 
     861""" 
     862    Z_MAX = 6.0    # maximum meaningful z-value 
     863    if z == 0.0: 
     864    x = 0.0 
     865    else: 
     866    y = 0.5 * math.fabs(z) 
     867    if y >= (Z_MAX*0.5): 
     868        x = 1.0 
     869    elif (y < 1.0): 
     870        w = y*y 
     871        x = ((((((((0.000124818987 * w 
     872            -0.001075204047) * w +0.005198775019) * w 
     873              -0.019198292004) * w +0.059054035642) * w 
     874            -0.151968751364) * w +0.319152932694) * w 
     875          -0.531923007300) * w +0.797884560593) * y * 2.0 
     876    else: 
     877        y = y - 2.0 
     878        x = (((((((((((((-0.000045255659 * y 
     879                 +0.000152529290) * y -0.000019538132) * y 
     880               -0.000676904986) * y +0.001390604284) * y 
     881             -0.000794620820) * y -0.002034254874) * y 
     882               +0.006549791214) * y -0.010557625006) * y 
     883             +0.011630447319) * y -0.009279453341) * y 
     884           +0.005353579108) * y -0.002141268741) * y 
     885         +0.000535310849) * y +0.999936657524 
     886    if z > 0.0: 
     887    prob = ((x+1.0)*0.5) 
     888    else: 
     889    prob = ((1.0-x)*0.5) 
     890    return prob 
     891 
     892    
  • orange/doc/Orange/rst/index.rst

    r7187 r7193  
    1515   orange.classify.svm 
    1616   orange.classify.tree 
     17   orange.classify.logreg 
    1718 
    1819Indices and tables 
  • orange/orngLR.py

    r6538 r7193  
    1 import orange 
    2 import orngCI 
    3 import math, os 
    4 from numpy import * 
    5 from numpy.linalg import * 
    6  
    7  
    8 ####################### 
    9 ## Print out methods ## 
    10 ####################### 
    11  
    12 def printOUT(classifier): 
    13     # print out class values 
    14     print 
    15     print "class attribute = " + classifier.domain.classVar.name 
    16     print "class values = " + str(classifier.domain.classVar.values) 
    17     print 
    18      
    19     # get the longest attribute name 
    20     longest=0 
    21     for at in classifier.continuizedDomain.attributes: 
    22         if len(at.name)>longest: 
    23             longest=len(at.name); 
    24  
    25     # print out the head 
    26     formatstr = "%"+str(longest)+"s %10s %10s %10s %10s %10s" 
    27     print formatstr % ("Attribute", "beta", "st. error", "wald Z", "P", "OR=exp(beta)") 
    28     print 
    29     formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f"     
    30     print formatstr % ("Intercept", classifier.beta[0], classifier.beta_se[0], classifier.wald_Z[0], classifier.P[0]) 
    31     formatstr = "%"+str(longest)+"s %10.2f %10.2f %10.2f %10.2f %10.2f"     
    32     for i in range(len(classifier.continuizedDomain.attributes)): 
    33         print formatstr % (classifier.continuizedDomain.attributes[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])) 
    34          
    35  
    36 def hasDiscreteValues(domain): 
    37     for at in domain.attributes: 
    38         if at.varType == orange.VarTypes.Discrete: 
    39             return 1 
    40     return 0 
    41  
    42 def LogRegLearner(examples = None, weightID=0, **kwds): 
    43     lr = LogRegLearnerClass(**kwds) 
    44     if examples: 
    45         return lr(examples, weightID) 
    46     else: 
    47         return lr 
    48  
    49 class LogRegLearnerClass(orange.Learner): 
    50     def __init__(self, removeSingular=0, fitter = None, **kwds): 
    51         self.__dict__.update(kwds) 
    52         self.removeSingular = removeSingular 
    53         self.fitter = None 
    54  
    55     def __call__(self, examples, weight=0): 
    56         imputer = getattr(self, "imputer", None) or None 
    57         if getattr(self, "removeMissing", 0): 
    58             examples = orange.Preprocessor_dropMissing(examples) 
    59 ##        if hasDiscreteValues(examples.domain): 
    60 ##            examples = createNoDiscTable(examples) 
    61         if not len(examples): 
    62             return None 
    63         if getattr(self, "stepwiseLR", 0): 
    64             addCrit = getattr(self, "addCrit", 0.2) 
    65             removeCrit = getattr(self, "removeCrit", 0.3) 
    66             numAttr = getattr(self, "numAttr", -1) 
    67             attributes = StepWiseFSS(examples, addCrit = addCrit, deleteCrit = removeCrit, imputer = imputer, numAttr = numAttr) 
    68             tmpDomain = orange.Domain(attributes, examples.domain.classVar) 
    69             tmpDomain.addmetas(examples.domain.getmetas()) 
    70             examples = examples.select(tmpDomain) 
    71         learner = orange.LogRegLearner() 
    72         learner.imputerConstructor = imputer 
    73         if imputer: 
    74             examples = self.imputer(examples)(examples) 
    75         examples = orange.Preprocessor_dropMissing(examples) 
    76         if self.fitter: 
    77             learner.fitter = self.fitter 
    78         if self.removeSingular: 
    79             lr = learner.fitModel(examples, weight) 
    80         else: 
    81             lr = learner(examples, weight) 
    82         while isinstance(lr, orange.Variable): 
    83             if isinstance(lr.getValueFrom, orange.ClassifierFromVar) and isinstance(lr.getValueFrom.transformer, orange.Discrete2Continuous): 
    84                 lr = lr.getValueFrom.variable 
    85             attributes = examples.domain.attributes[:] 
    86             if lr in attributes: 
    87                 attributes.remove(lr) 
    88             else: 
    89                 attributes.remove(lr.getValueFrom.variable) 
    90             newDomain = orange.Domain(attributes, examples.domain.classVar) 
    91             newDomain.addmetas(examples.domain.getmetas()) 
    92             examples = examples.select(newDomain) 
    93             lr = learner.fitModel(examples, weight) 
    94         return lr 
    95  
    96  
    97  
    98 def Univariate_LogRegLearner(examples=None, **kwds): 
    99     learner = apply(Univariate_LogRegLearner_Class, (), kwds) 
    100     if examples: 
    101         return learner(examples) 
    102     else: 
    103         return learner 
    104  
    105 class Univariate_LogRegLearner_Class(orange.Learner): 
    106     def __init__(self, **kwds): 
    107         self.__dict__.update(kwds) 
    108  
    109     def __call__(self, examples): 
    110         examples = createFullNoDiscTable(examples) 
    111         classifiers = map(lambda x: LogRegLearner(orange.Preprocessor_dropMissing(examples.select(orange.Domain(x, examples.domain.classVar)))), examples.domain.attributes) 
    112         maj_classifier = LogRegLearner(orange.Preprocessor_dropMissing(examples.select(orange.Domain(examples.domain.classVar)))) 
    113         beta = [maj_classifier.beta[0]] + [x.beta[1] for x in classifiers] 
    114         beta_se = [maj_classifier.beta_se[0]] + [x.beta_se[1] for x in classifiers] 
    115         P = [maj_classifier.P[0]] + [x.P[1] for x in classifiers] 
    116         wald_Z = [maj_classifier.wald_Z[0]] + [x.wald_Z[1] for x in classifiers] 
    117         domain = examples.domain 
    118  
    119         return Univariate_LogRegClassifier(beta = beta, beta_se = beta_se, P = P, wald_Z = wald_Z, domain = domain) 
    120  
    121 class Univariate_LogRegClassifier(orange.Classifier): 
    122     def __init__(self, **kwds): 
    123         self.__dict__.update(kwds) 
    124  
    125     def __call__(self, example, resultType = orange.GetValue): 
    126         # classification not implemented yet. For now its use is only to provide regression coefficients and its statistics 
    127         pass 
    128      
    129  
    130 def LogRegLearner_getPriors(examples = None, weightID=0, **kwds): 
    131     lr = LogRegLearnerClass_getPriors(**kwds) 
    132     if examples: 
    133         return lr(examples, weightID) 
    134     else: 
    135         return lr 
    136  
    137 class LogRegLearnerClass_getPriors(object): 
    138     def __init__(self, removeSingular=0, **kwds): 
    139         self.__dict__.update(kwds) 
    140         self.removeSingular = removeSingular 
    141     def __call__(self, examples, weight=0): 
    142         # next function changes data set to a extended with unknown values  
    143         def createLogRegExampleTable(data, weightID): 
    144             setsOfData = [] 
    145             for at in data.domain.attributes: 
    146                 # za vsak atribut kreiraj nov newExampleTable newData 
    147                 # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable 
    148                 if at.varType == orange.VarTypes.Continuous: 
    149                     atDisc = orange.FloatVariable(at.name + "Disc") 
    150                     newDomain = orange.Domain(data.domain.attributes+[atDisc,data.domain.classVar]) 
    151                     newDomain.addmetas(data.domain.getmetas()) 
    152                     newData = orange.ExampleTable(newDomain,data) 
    153                     altData = orange.ExampleTable(newDomain,data) 
    154                     for i,d in enumerate(newData): 
    155                         d[atDisc] = 0 
    156                         d[weightID] = 1*data[i][weightID] 
    157                     for i,d in enumerate(altData): 
    158                         d[atDisc] = 1 
    159                         d[at] = 0 
    160                         d[weightID] = 0.000001*data[i][weightID] 
    161                 elif at.varType == orange.VarTypes.Discrete: 
    162                 # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X" 
    163                     atNew = orange.EnumVariable(at.name, values = at.values + [at.name+"X"]) 
    164                     newDomain = orange.Domain(filter(lambda x: x!=at, data.domain.attributes)+[atNew,data.domain.classVar]) 
    165                     newDomain.addmetas(data.domain.getmetas()) 
    166                     newData = orange.ExampleTable(newDomain,data) 
    167                     altData = orange.ExampleTable(newDomain,data) 
    168                     for i,d in enumerate(newData): 
    169                         d[atNew] = data[i][at] 
    170                         d[weightID] = 1*data[i][weightID] 
    171                     for i,d in enumerate(altData): 
    172                         d[atNew] = at.name+"X" 
    173                         d[weightID] = 0.000001*data[i][weightID] 
    174                 newData.extend(altData) 
    175                 setsOfData.append(newData) 
    176             return setsOfData 
    177                    
    178         learner = LogRegLearner(imputer = orange.ImputerConstructor_average(), removeSingular = self.removeSingular) 
    179         # get Original Model 
    180         orig_model = learner(examples,weight) 
    181         if orig_model.fit_status: 
    182             print "Warning: model did not converge" 
    183  
    184         # get extended Model (you should not change data) 
    185         if weight == 0: 
    186             weight = orange.newmetaid() 
    187             examples.addMetaAttribute(weight, 1.0) 
    188         extended_set_of_examples = createLogRegExampleTable(examples, weight) 
    189         extended_models = [learner(extended_examples, weight) \ 
    190                            for extended_examples in extended_set_of_examples] 
    191  
    192 ##        print examples[0] 
    193 ##        printOUT(orig_model) 
    194 ##        print orig_model.domain 
    195 ##        print orig_model.beta 
    196 ##        print orig_model.beta[orig_model.continuizedDomain.attributes[-1]] 
    197 ##        for i,m in enumerate(extended_models): 
    198 ##            print examples.domain.attributes[i] 
    199 ##            printOUT(m) 
    200              
    201          
    202         # izracunas odstopanja 
    203         # get sum of all betas 
    204         beta = 0 
    205         betas_ap = [] 
    206         for m in extended_models: 
    207             beta_add = m.beta[m.continuizedDomain.attributes[-1]] 
    208             betas_ap.append(beta_add) 
    209             beta = beta + beta_add 
    210          
    211         # substract it from intercept 
    212         #print "beta", beta 
    213         logistic_prior = orig_model.beta[0]+beta 
    214          
    215         # compare it to bayes prior 
    216         bayes = orange.BayesLearner(examples) 
    217         bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0]) 
    218  
    219         # normalize errors 
    220 ##        print "bayes", bayes_prior 
    221 ##        print "lr", orig_model.beta[0] 
    222 ##        print "lr2", logistic_prior 
    223 ##        print "dist", orange.Distribution(examples.domain.classVar,examples) 
    224 ##        print "prej", betas_ap 
    225  
    226         # error normalization - to avoid errors due to assumption of independence of unknown values 
    227         dif = bayes_prior - logistic_prior 
    228         positives = sum(filter(lambda x: x>=0, betas_ap)) 
    229         negatives = -sum(filter(lambda x: x<0, betas_ap)) 
    230         if not negatives == 0: 
    231             kPN = positives/negatives 
    232             diffNegatives = dif/(1+kPN) 
    233             diffPositives = kPN*diffNegatives 
    234             kNegatives = (negatives-diffNegatives)/negatives 
    235             kPositives = positives/(positives-diffPositives) 
    236     ##        print kNegatives 
    237     ##        print kPositives 
    238  
    239             for i,b in enumerate(betas_ap): 
    240                 if b<0: betas_ap[i]*=kNegatives 
    241                 else: betas_ap[i]*=kPositives 
    242         #print "potem", betas_ap 
    243  
    244         # vrni originalni model in pripadajoce apriorne niclele 
    245         return (orig_model, betas_ap) 
    246         #return (bayes_prior,orig_model.beta[examples.domain.classVar],logistic_prior) 
    247  
    248 class LogRegLearnerClass_getPriors_OneTable: 
    249     def __init__(self, removeSingular=0, **kwds): 
    250         self.__dict__.update(kwds) 
    251         self.removeSingular = removeSingular 
    252     def __call__(self, examples, weight=0): 
    253         # next function changes data set to a extended with unknown values  
    254         def createLogRegExampleTable(data, weightID): 
    255             finalData = orange.ExampleTable(data) 
    256             origData = orange.ExampleTable(data) 
    257             for at in data.domain.attributes: 
    258                 # za vsak atribut kreiraj nov newExampleTable newData 
    259                 # v dataOrig, dataFinal in newData dodaj nov atribut -- continuous variable 
    260                 if at.varType == orange.VarTypes.Continuous: 
    261                     atDisc = orange.FloatVariable(at.name + "Disc") 
    262                     newDomain = orange.Domain(origData.domain.attributes+[atDisc,data.domain.classVar]) 
    263                     newDomain.addmetas(newData.domain.getmetas()) 
    264                     finalData = orange.ExampleTable(newDomain,finalData) 
    265                     newData = orange.ExampleTable(newDomain,origData) 
    266                     origData = orange.ExampleTable(newDomain,origData) 
    267                     for d in origData: 
    268                         d[atDisc] = 0 
    269                     for d in finalData: 
    270                         d[atDisc] = 0 
    271                     for i,d in enumerate(newData): 
    272                         d[atDisc] = 1 
    273                         d[at] = 0 
    274                         d[weightID] = 100*data[i][weightID] 
    275                          
    276                 elif at.varType == orange.VarTypes.Discrete: 
    277                 # v dataOrig, dataFinal in newData atributu "at" dodaj ee  eno  vreednost, ki ima vrednost kar  ime atributa +  "X" 
    278                     atNew = orange.EnumVariable(at.name, values = at.values + [at.name+"X"]) 
    279                     newDomain = orange.Domain(filter(lambda x: x!=at, origData.domain.attributes)+[atNew,origData.domain.classVar]) 
    280                     newDomain.addmetas(origData.domain.getmetas()) 
    281                     temp_finalData = orange.ExampleTable(finalData) 
    282                     finalData = orange.ExampleTable(newDomain,finalData) 
    283                     newData = orange.ExampleTable(newDomain,origData) 
    284                     temp_origData = orange.ExampleTable(origData) 
    285                     origData = orange.ExampleTable(newDomain,origData) 
    286                     for i,d in enumerate(origData): 
    287                         d[atNew] = temp_origData[i][at] 
    288                     for i,d in enumerate(finalData): 
    289                         d[atNew] = temp_finalData[i][at]                         
    290                     for i,d in enumerate(newData): 
    291                         d[atNew] = at.name+"X" 
    292                         d[weightID] = 10*data[i][weightID] 
    293                 finalData.extend(newData) 
    294             return finalData 
    295                    
    296         learner = LogRegLearner(imputer = orange.ImputerConstructor_average(), removeSingular = self.removeSingular) 
    297         # get Original Model 
    298         orig_model = learner(examples,weight) 
    299  
    300         # get extended Model (you should not change data) 
    301         if weight == 0: 
    302             weight = orange.newmetaid() 
    303             examples.addMetaAttribute(weight, 1.0) 
    304         extended_examples = createLogRegExampleTable(examples, weight) 
    305         extended_model = learner(extended_examples, weight) 
    306  
    307 ##        print examples[0] 
    308 ##        printOUT(orig_model) 
    309 ##        print orig_model.domain 
    310 ##        print orig_model.beta 
    311  
    312 ##        printOUT(extended_model)         
    313         # izracunas odstopanja 
    314         # get sum of all betas 
    315         beta = 0 
    316         betas_ap = [] 
    317         for m in extended_models: 
    318             beta_add = m.beta[m.continuizedDomain.attributes[-1]] 
    319             betas_ap.append(beta_add) 
    320             beta = beta + beta_add 
    321          
    322         # substract it from intercept 
    323         #print "beta", beta 
    324         logistic_prior = orig_model.beta[0]+beta 
    325          
    326         # compare it to bayes prior 
    327         bayes = orange.BayesLearner(examples) 
    328         bayes_prior = math.log(bayes.distribution[1]/bayes.distribution[0]) 
    329  
    330         # normalize errors 
    331         #print "bayes", bayes_prior 
    332         #print "lr", orig_model.beta[0] 
    333         #print "lr2", logistic_prior 
    334         #print "dist", orange.Distribution(examples.domain.classVar,examples) 
    335         k = (bayes_prior-orig_model.beta[0])/(logistic_prior-orig_model.beta[0]) 
    336         #print "prej", betas_ap 
    337         betas_ap = [k*x for x in betas_ap]                 
    338         #print "potem", betas_ap 
    339  
    340         # vrni originalni model in pripadajoce apriorne niclele 
    341         return (orig_model, betas_ap) 
    342         #return (bayes_prior,orig_model.beta[data.domain.classVar],logistic_prior) 
    343  
    344  
    345 ###################################### 
    346 #### Fitters for logistic regression (logreg) learner #### 
    347 ###################################### 
    348  
    349 def Pr(x, betas): 
    350     k = math.exp(dot(x, betas)) 
    351     return k / (1+k) 
    352  
    353 def lh(x,y,betas): 
    354     llh = 0.0 
    355     for i,x_i in enumerate(x): 
    356         pr = Pr(x_i,betas) 
    357         llh += y[i]*log(max(pr,1e-6)) + (1-y[i])*log(max(1-pr,1e-6)) 
    358     return llh 
    359  
    360  
    361 def diag(vector): 
    362     mat = identity(len(vector), Float) 
    363     for i,v in enumerate(vector): 
    364         mat[i][i] = v 
    365     return mat 
    366      
    367 class simpleFitter(orange.LogRegFitter): 
    368     def __init__(self, penalty=0, se_penalty = False): 
    369         self.penalty = penalty 
    370         self.se_penalty = se_penalty 
    371     def __call__(self, data, weight=0): 
    372         ml = data.native(0) 
    373         for i in range(len(data.domain.attributes)): 
    374           a = data.domain.attributes[i] 
    375           if a.varType == orange.VarTypes.Discrete: 
    376             for m in ml: 
    377               m[i] = a.values.index(m[i]) 
    378         for m in ml: 
    379           m[-1] = data.domain.classVar.values.index(m[-1]) 
    380         Xtmp = array(ml) 
    381         y = Xtmp[:,-1]   # true probabilities (1's or 0's) 
    382         one = reshape(array([1]*len(data)), (len(data),1)) # intercept column 
    383         X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data 
    384  
    385         betas = array([0.0] * (len(data.domain.attributes)+1)) 
    386         oldBetas = array([1.0] * (len(data.domain.attributes)+1)) 
    387         N = len(data) 
    388  
    389         pen_matrix = array([self.penalty] * (len(data.domain.attributes)+1)) 
    390         if self.se_penalty: 
    391             p = array([Pr(X[i], betas) for i in range(len(data))]) 
    392             W = identity(len(data), Float) 
    393             pp = p * (1.0-p) 
    394             for i in range(N): 
    395                 W[i,i] = pp[i] 
    396             se = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))))) 
    397             for i,p in enumerate(pen_matrix): 
    398                 pen_matrix[i] *= se[i] 
    399         # predict the probability for an instance, x and betas are vectors 
    400         # start the computation 
    401         likelihood = 0. 
    402         likelihood_new = 1. 
    403         while abs(likelihood - likelihood_new)>1e-5: 
    404             likelihood = likelihood_new 
    405             oldBetas = betas 
    406             p = array([Pr(X[i], betas) for i in range(len(data))]) 
    407  
    408             W = identity(len(data), Float) 
    409             pp = p * (1.0-p) 
    410             for i in range(N): 
    411                 W[i,i] = pp[i] 
    412  
    413             WI = inverse(W) 
    414             z = matrixmultiply(X, betas) + matrixmultiply(WI, y - p) 
    415  
    416             tmpA = inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))+diag(pen_matrix)) 
    417             tmpB = matrixmultiply(transpose(X), y-p) 
    418             betas = oldBetas + matrixmultiply(tmpA,tmpB) 
    419 #            betaTemp = matrixmultiply(matrixmultiply(matrixmultiply(matrixmultiply(tmpA,transpose(X)),W),X),oldBetas) 
    420 #            print betaTemp 
    421 #            tmpB = matrixmultiply(transpose(X), matrixmultiply(W, z)) 
    422 #            betas = matrixmultiply(tmpA, tmpB) 
    423             likelihood_new = lh(X,y,betas)-self.penalty*sum([b*b for b in betas]) 
    424             print likelihood_new 
    425  
    426              
    427              
    428 ##        XX = sqrt(diagonal(inverse(matrixmultiply(transpose(X),X)))) 
    429 ##        yhat = array([Pr(X[i], betas) for i in range(len(data))]) 
    430 ##        ss = sum((y - yhat) ** 2) / (N - len(data.domain.attributes) - 1) 
    431 ##        sigma = math.sqrt(ss) 
    432         p = array([Pr(X[i], betas) for i in range(len(data))]) 
    433         W = identity(len(data), Float) 
    434         pp = p * (1.0-p) 
    435         for i in range(N): 
    436             W[i,i] = pp[i] 
    437         diXWX = sqrt(diagonal(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))))) 
    438         xTemp = matrixmultiply(matrixmultiply(inverse(matrixmultiply(transpose(X), matrixmultiply(W, X))),transpose(X)),y) 
    439         beta = [] 
    440         beta_se = [] 
    441         print "likelihood ridge", likelihood 
    442         for i in range(len(betas)): 
    443             beta.append(betas[i]) 
    444             beta_se.append(diXWX[i]) 
    445         return (self.OK, beta, beta_se, 0) 
    446  
    447 def Pr_bx(bx): 
    448     if bx > 35: 
    449         return 1 
    450     if bx < -35: 
    451         return 0 
    452     return exp(bx)/(1+exp(bx)) 
    453  
    454 class bayesianFitter(orange.LogRegFitter): 
    455     def __init__(self, penalty=0, anch_examples=[], tau = 0): 
    456         self.penalty = penalty 
    457         self.anch_examples = anch_examples 
    458         self.tau = tau 
    459  
    460     def createArrayData(self,data): 
    461         if not len(data): 
    462             return (array([]),array([])) 
    463         # convert data to numeric 
    464         ml = data.native(0) 
    465         for i,a in enumerate(data.domain.attributes): 
    466           if a.varType == orange.VarTypes.Discrete: 
    467             for m in ml: 
    468               m[i] = a.values.index(m[i]) 
    469         for m in ml: 
    470           m[-1] = data.domain.classVar.values.index(m[-1]) 
    471         Xtmp = array(ml) 
    472         y = Xtmp[:,-1]   # true probabilities (1's or 0's) 
    473         one = reshape(array([1]*len(data)), (len(data),1)) # intercept column 
    474         X=concatenate((one, Xtmp[:,:-1]),1)  # intercept first, then data 
    475         return (X,y) 
    476      
    477     def __call__(self, data, weight=0): 
    478         (X,y)=self.createArrayData(data) 
    479  
    480         exTable = orange.ExampleTable(data.domain) 
    481         for id,ex in self.anch_examples: 
    482             exTable.extend(orange.ExampleTable(ex,data.domain)) 
    483         (X_anch,y_anch)=self.createArrayData(exTable) 
    484  
    485         betas = array([0.0] * (len(data.domain.attributes)+1)) 
    486  
    487         likelihood,betas = self.estimateBeta(X,y,betas,[0]*(len(betas)),X_anch,y_anch) 
    488  
    489         # get attribute groups atGroup = [(startIndex, number of values), ...) 
    490         ats = data.domain.attributes 
    491         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:] 
    492         atGroup=[[0,0]] 
    493         for v_i,v in enumerate(atVec): 
    494             if v[1]==0: atGroup[-1][1]+=1 
    495             else:       atGroup.append([v_i,1]) 
    496          
    497         # compute zero values for attributes 
    498         sumB = 0. 
    499         for ag in atGroup: 
    500             X_temp = concatenate((X[:,:ag[0]+1],X[:,ag[0]+1+ag[1]:]),1) 
    501             if X_anch: 
    502                 X_anch_temp = concatenate((X_anch[:,:ag[0]+1],X_anch[:,ag[0]+1+ag[1]:]),1) 
    503             else: X_anch_temp = X_anch 
    504 ##            print "1", concatenate((betas[:i+1],betas[i+2:])) 
    505 ##            print "2", betas 
    506             likelihood_temp,betas_temp=self.estimateBeta(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) 
    507             print "finBetas", betas, betas_temp 
    508             print "betas", betas[0], betas_temp[0] 
    509             sumB += betas[0]-betas_temp[0] 
    510         apriori = orange.Distribution(data.domain.classVar, data) 
    511         aprioriProb = apriori[0]/apriori.abs 
    512          
    513         print "koncni rezultat", sumB, math.log((1-aprioriProb)/aprioriProb), betas[0] 
    514              
    515         beta = [] 
    516         beta_se = [] 
    517         print "likelihood2", likelihood 
    518         for i in range(len(betas)): 
    519             beta.append(betas[i]) 
    520             beta_se.append(0.0) 
    521         return (self.OK, beta, beta_se, 0) 
    522  
    523       
    524          
    525     def estimateBeta(self,X,y,betas,const_betas,X_anch,y_anch): 
    526         N,N_anch = len(y),len(y_anch) 
    527         r,r_anch = array([dot(X[i], betas) for i in range(N)]),\ 
    528                    array([dot(X_anch[i], betas) for i in range(N_anch)]) 
    529         p    = array([Pr_bx(ri) for ri in r]) 
    530         X_sq = X*X 
    531  
    532         max_delta      = [1.]*len(const_betas) 
    533         likelihood     = -1.e+10 
    534         likelihood_new = -1.e+9 
    535         while abs(likelihood - likelihood_new)>0.01 and max(max_delta)>0.01: 
    536             likelihood = likelihood_new 
    537             print likelihood 
    538             betas_temp = [b for b in betas] 
    539             for j in range(len(betas)): 
    540                 if const_betas[j]: continue 
    541                 dl = matrixmultiply(X[:,j],transpose(y-p)) 
    542                 for xi,x in enumerate(X_anch): 
    543                     dl += self.penalty*x[j]*(y_anch[xi] - Pr_bx(r_anch[xi]*self.penalty)) 
    544  
    545                 ddl = matrixmultiply(X_sq[:,j],transpose(p*(1-p))) 
    546                 for xi,x in enumerate(X_anch): 
    547                     ddl += self.penalty*x[j]*Pr_bx(r[xi]*self.penalty)*(1-Pr_bx(r[xi]*self.penalty)) 
    548  
    549                 if j==0: 
    550                     dv = dl/max(ddl,1e-6) 
    551                 elif betas[j] == 0: # special handling due to non-defined first and second derivatives 
    552                     dv = (dl-self.tau)/max(ddl,1e-6) 
    553                     if dv < 0: 
    554                         dv = (dl+self.tau)/max(ddl,1e-6) 
    555                         if dv > 0: 
    556                             dv = 0 
    557                 else: 
    558                     dl -= sign(betas[j])*self.tau 
    559                     dv = dl/max(ddl,1e-6) 
    560                     if not sign(betas[j] + dv) == sign(betas[j]): 
    561                         dv = -betas[j] 
    562                 dv = min(max(dv,-max_delta[j]),max_delta[j]) 
    563                 r+= X[:,j]*dv 
    564                 p = array([Pr_bx(ri) for ri in r]) 
    565                 if N_anch: 
    566                     r_anch+=X_anch[:,j]*dv 
    567                 betas[j] += dv 
    568                 max_delta[j] = max(2*abs(dv),max_delta[j]/2) 
    569             likelihood_new = lh(X,y,betas) 
    570             for xi,x in enumerate(X_anch): 
    571                 try: 
    572                     likelihood_new += y_anch[xi]*r_anch[xi]*self.penalty-log(1+exp(r_anch[xi]*self.penalty)) 
    573                 except: 
    574                     likelihood_new += r_anch[xi]*self.penalty*(y_anch[xi]-1) 
    575             likelihood_new -= sum([abs(b) for b in betas[1:]])*self.tau 
    576             if likelihood_new < likelihood: 
    577                 max_delta = [md/4 for md in max_delta] 
    578                 likelihood_new = likelihood 
    579                 likelihood = likelihood_new + 1. 
    580                 betas = [b for b in betas_temp] 
    581         print "betas", betas 
    582         print "init_like", likelihood_new 
    583         print "pure_like", lh(X,y,betas) 
    584         return (likelihood,betas) 
    585      
    586 ############################################################ 
    587 ####  Feature subset selection for logistic regression  #### 
    588 ############################################################ 
    589  
    590  
    591 def StepWiseFSS(examples = None, **kwds): 
    592     """ 
    593       Constructs and returns a new set of examples that includes a 
    594       class and attributes selected by stepwise logistic regression. This is an 
    595       implementation of algorithm described in [Hosmer and Lemeshow, Applied Logistic Regression, 2000] 
    596  
    597       examples: data set (ExampleTable)      
    598       addCrit: "Alpha" level to judge if variable has enough importance to be added in the new set. (e.g. if addCrit is 0.2, then attribute is added if its P is lower than 0.2) 
    599       deleteCrit: Similar to addCrit, just that it is used at backward elimination. It should be higher than addCrit! 
    600       numAttr: maximum number of selected attributes, use -1 for infinity 
    601     """ 
    602  
    603     fss = apply(StepWiseFSS_class, (), kwds) 
    604     if examples is not None: 
    605         return fss(examples) 
    606     else: 
    607         return fss 
    608  
    609 def getLikelihood(fitter, examples): 
    610     res = fitter(examples) 
    611     if res[0] in [fitter.OK]: #, fitter.Infinity, fitter.Divergence]: 
    612        status, beta, beta_se, likelihood = res 
    613        if sum([abs(b) for b in beta])<sum([abs(b) for b in beta_se]): 
    614            return -100*len(examples) 
    615        return likelihood 
    616     else: 
    617        return -100*len(examples) 
    618          
    619  
    620  
    621 class StepWiseFSS_class(orange.Learner): 
    622   def __init__(self, addCrit=0.2, deleteCrit=0.3, numAttr = -1, **kwds): 
    623     self.__dict__.update(kwds) 
    624     self.addCrit = addCrit 
    625     self.deleteCrit = deleteCrit 
    626     self.numAttr = numAttr 
    627   def __call__(self, examples): 
    628     if getattr(self, "imputer", 0): 
    629         examples = self.imputer(examples)(examples) 
    630     if getattr(self, "removeMissing", 0): 
    631         examples = orange.Preprocessor_dropMissing(examples) 
    632     continuizer = orange.DomainContinuizer(zeroBased=1,continuousTreatment=orange.DomainContinuizer.Leave, 
    633                                            multinomialTreatment = orange.DomainContinuizer.FrequentIsBase, 
    634                                            classTreatment = orange.DomainContinuizer.Ignore) 
    635     attr = [] 
    636     remain_attr = examples.domain.attributes[:] 
    637  
    638     # get LL for Majority Learner  
    639     tempDomain = orange.Domain(attr,examples.domain.classVar) 
    640     #tempData  = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
    641     tempData  = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
    642  
    643     ll_Old = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
    644     ll_Best = -1000000 
    645     length_Old = float(len(tempData)) 
    646  
    647     stop = 0 
    648     while not stop: 
    649         # LOOP until all variables are added or no further deletion nor addition of attribute is possible 
    650         worstAt = None 
    651         # if there are more than 1 attribute then perform backward elimination 
    652         if len(attr) >= 2: 
    653             minG = 1000 
    654             worstAt = attr[0] 
    655             ll_Best = ll_Old 
    656             length_Best = length_Old 
    657             for at in attr: 
    658                 # check all attribute whether its presence enough increases LL? 
    659  
    660                 tempAttr = filter(lambda x: x!=at, attr) 
    661                 tempDomain = orange.Domain(tempAttr,examples.domain.classVar) 
    662                 tempDomain.addmetas(examples.domain.getmetas()) 
    663                 # domain, calculate P for LL improvement. 
    664                 tempDomain  = continuizer(orange.Preprocessor_dropMissing(examples.select(tempDomain))) 
    665                 tempData = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
    666  
    667                 ll_Delete = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
    668                 length_Delete = float(len(tempData)) 
    669                 length_Avg = (length_Delete + length_Old)/2.0 
    670  
    671                 G=-2*length_Avg*(ll_Delete/length_Delete-ll_Old/length_Old) 
    672  
    673                 # set new worst attribute                 
    674                 if G<minG: 
    675                     worstAt = at 
    676                     minG=G 
    677                     ll_Best = ll_Delete 
    678                     length_Best = length_Delete 
    679             # deletion of attribute 
    680              
    681             if worstAt.varType==orange.VarTypes.Continuous: 
    682                 P=lchisqprob(minG,1); 
    683             else: 
    684                 P=lchisqprob(minG,len(worstAt.values)-1); 
    685             if P>=self.deleteCrit: 
    686                 attr.remove(worstAt) 
    687                 remain_attr.append(worstAt) 
    688                 nodeletion=0 
    689                 ll_Old = ll_Best 
    690                 length_Old = length_Best 
    691             else: 
    692                 nodeletion=1 
    693         else: 
    694             nodeletion = 1 
    695             # END OF DELETION PART 
    696              
    697         # if enough attributes has been chosen, stop the procedure 
    698         if self.numAttr>-1 and len(attr)>=self.numAttr: 
    699             remain_attr=[] 
    700           
    701         # for each attribute in the remaining 
    702         maxG=-1 
    703         ll_Best = ll_Old 
    704         length_Best = length_Old 
    705         bestAt = None 
    706         for at in remain_attr: 
    707             tempAttr = attr + [at] 
    708             tempDomain = orange.Domain(tempAttr,examples.domain.classVar) 
    709             tempDomain.addmetas(examples.domain.getmetas()) 
    710             # domain, calculate P for LL improvement. 
    711             tempDomain  = continuizer(orange.Preprocessor_dropMissing(examples.select(tempDomain))) 
    712             tempData = orange.Preprocessor_dropMissing(examples.select(tempDomain)) 
    713             ll_New = getLikelihood(orange.LogRegFitter_Cholesky(), tempData) 
    714  
    715             length_New = float(len(tempData)) # get number of examples in tempData to normalize likelihood 
    716  
    717             # P=PR(CHI^2>G), G=-2(L(0)-L(1))=2(E(0)-E(1)) 
    718             length_avg = (length_New + length_Old)/2 
    719             G=-2*length_avg*(ll_Old/length_Old-ll_New/length_New); 
    720             if G>maxG: 
    721                 bestAt = at 
    722                 maxG=G 
    723                 ll_Best = ll_New 
    724                 length_Best = length_New 
    725         if not bestAt: 
    726             stop = 1 
    727             continue 
    728          
    729         if bestAt.varType==orange.VarTypes.Continuous: 
    730             P=lchisqprob(maxG,1); 
    731         else: 
    732             P=lchisqprob(maxG,len(bestAt.values)-1); 
    733         # Add attribute with smallest P to attributes(attr) 
    734         if P<=self.addCrit: 
    735             attr.append(bestAt) 
    736             remain_attr.remove(bestAt) 
    737             ll_Old = ll_Best 
    738             length_Old = length_Best 
    739  
    740         if (P>self.addCrit and nodeletion) or (bestAt == worstAt): 
    741             stop = 1 
    742  
    743     return attr 
    744  
    745  
    746 def StepWiseFSS_Filter(examples = None, **kwds): 
    747     """ 
    748         check function StepWiseFSS() 
    749     """ 
    750  
    751     filter = apply(StepWiseFSS_Filter_class, (), kwds) 
    752     if examples is not None: 
    753         return filter(examples) 
    754     else: 
    755         return filter 
    756  
    757  
    758 class StepWiseFSS_Filter_class(object): 
    759     def __init__(self, addCrit=0.2, deleteCrit=0.3, numAttr = -1): 
    760         self.addCrit = addCrit 
    761         self.deleteCrit = deleteCrit 
    762         self.numAttr = numAttr 
    763     def __call__(self, examples): 
    764         attr = StepWiseFSS(examples, addCrit=self.addCrit, deleteCrit = self.deleteCrit, numAttr = self.numAttr) 
    765         return examples.select(orange.Domain(attr, examples.domain.classVar)) 
    766                  
    767  
    768 #################################### 
    769 ####  PROBABILITY CALCULATIONS  #### 
    770 #################################### 
    771  
    772 def lchisqprob(chisq,df): 
    773     """ 
    774 Returns the (1-tailed) probability value associated with the provided 
    775 chi-square value and df.  Adapted from chisq.c in Gary Perlman's |Stat. 
    776  
    777 Usage:   lchisqprob(chisq,df) 
    778 """ 
    779     BIG = 20.0 
    780     def ex(x): 
    781         BIG = 20.0 
    782         if x < -BIG: 
    783             return 0.0 
    784         else: 
    785             return math.exp(x) 
    786     if chisq <=0 or df < 1: 
    787         return 1.0 
    788     a = 0.5 * chisq 
    789     if df%2 == 0: 
    790         even = 1 
    791     else: 
    792         even = 0 
    793     if df > 1: 
    794         y = ex(-a) 
    795     if even: 
    796         s = y 
    797     else: 
    798         s = 2.0 * zprob(-math.sqrt(chisq)) 
    799     if (df > 2): 
    800         chisq = 0.5 * (df - 1.0) 
    801         if even: 
    802             z = 1.0 
    803         else: 
    804             z = 0.5 
    805         if a > BIG: 
    806             if even: 
    807                 e = 0.0 
    808             else: 
    809                 e = math.log(math.sqrt(math.pi)) 
    810             c = math.log(a) 
    811             while (z <= chisq): 
    812                 e = math.log(z) + e 
    813                 s = s + ex(c*z-a-e) 
    814                 z = z + 1.0 
    815             return s 
    816         else: 
    817             if even: 
    818                 e = 1.0 
    819             else: 
    820                 e = 1.0 / math.sqrt(math.pi) / math.sqrt(a) 
    821             c = 0.0 
    822             while (z <= chisq): 
    823                 e = e * (a/float(z)) 
    824                 c = c + e 
    825                 z = z + 1.0 
    826             return (c*y+s) 
    827     else: 
    828         return s 
    829  
    830  
    831 def zprob(z): 
    832     """ 
    833 Returns the area under the normal curve 'to the left of' the given z value. 
    834 Thus,  
    835     for z<0, zprob(z) = 1-tail probability 
    836     for z>0, 1.0-zprob(z) = 1-tail probability 
    837     for any z, 2.0*(1.0-zprob(abs(z))) = 2-tail probability 
    838 Adapted from z.c in Gary Perlman's |Stat. 
    839  
    840 Usage:   lzprob(z) 
    841 """ 
    842     Z_MAX = 6.0    # maximum meaningful z-value 
    843     if z == 0.0: 
    844     x = 0.0 
    845     else: 
    846     y = 0.5 * math.fabs(z) 
    847     if y >= (Z_MAX*0.5): 
    848         x = 1.0 
    849     elif (y < 1.0): 
    850         w = y*y 
    851         x = ((((((((0.000124818987 * w 
    852             -0.001075204047) * w +0.005198775019) * w 
    853               -0.019198292004) * w +0.059054035642) * w 
    854             -0.151968751364) * w +0.319152932694) * w 
    855           -0.531923007300) * w +0.797884560593) * y * 2.0 
    856     else: 
    857         y = y - 2.0 
    858         x = (((((((((((((-0.000045255659 * y 
    859                  +0.000152529290) * y -0.000019538132) * y 
    860                -0.000676904986) * y +0.001390604284) * y 
    861              -0.000794620820) * y -0.002034254874) * y 
    862                +0.006549791214) * y -0.010557625006) * y 
    863              +0.011630447319) * y -0.009279453341) * y 
    864            +0.005353579108) * y -0.002141268741) * y 
    865          +0.000535310849) * y +0.999936657524 
    866     if z > 0.0: 
    867     prob = ((x+1.0)*0.5) 
    868     else: 
    869     prob = ((1.0-x)*0.5) 
    870     return prob 
    871  
    872     
     1from Orange.classify.logreg import * 
Note: See TracChangeset for help on using the changeset viewer.