Changeset 7206:f4153e99e376 in orange


Ignore:
Timestamp:
02/02/11 17:32:59 (3 years ago)
Author:
jzbontar <jure.zbontar@…>
Branch:
default
Convert:
b6ec0fb1d9c0f847a5b89f52cf8e5c36876ea0d4
Message:
 
File:
1 edited

Legend:

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

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