CS535D Project: Bayesian Logistic Regression through Auxiliary Variables

Mark Schmidt

Description

Niveau: Supérieur, Licence, Bac+1
CS535D Project: Bayesian Logistic Regression through Auxiliary Variables Mark Schmidt Abstract This project deals with the estimation of Logistic Regression parameters. We first review the binary logistic regression model and the multinomial extension, including standard MAP parameter estimation with a Gaussian prior. We then turn to the case of Bayesian Logistic Regression under this same prior. We review the cannonical approach of performing Bayesian Probit Regression through auxiliary variables, and extensions of this technique to Bayesian Logistic Regression and Bayesian Multinomial Regression. We then turn to the task of feature selection, outlining a trans-dimensional MCMC approach to variable selection in Bayesian Logistic Regression. Finally, we turn to the case of estimating MAP parameters and performing Bayesian Logistic Regression under L1 penalties and other sparsity promoting priors. 1 Introduction In this project, we examined the highly popular Logistic Regression model. This model has tradition- ally been appealing due to its performance in classification, the potential to use its outputs as probabilitic estimates since they are in the range [0, 1], and the interpretation of the coefficients in terms of the 'log- odds' ratio [1]. It is especially popular in biostatistical applications where binary classification tasks occur frequently [1]. In this first part of the report, we review this model, its multi-class generalization, and standard methods of performing maximum likelihood (ML) or maximum a posteriori (MAP) para- meter estimation under a zero-mean Gaussian prior for the regression coefficients.

• gaussian


• parameter estimation


• multinomial


• auxiliary variable


• bayesian logistic


• regression


• binary probit


• function

Sujets

• Schmidt
• Gaussian
• Estimation theory
• Multinomial
• Régression
• Function

Informations

• Publié par: profil-ontoe-2012
• Nombre de lectures: 18
• Langue: English

Extrait

CS535DProject:BayesianLogisticRegressionthroughAuxiliary Variables MarkSchmidt Abstract ThisprojectdealswiththeestimationofLogisticRegressionparameters.Werstreviewthebinary logisticregressionmodelandthemultinomialextension,includingstandardMAPparameterestimation withaGaussianprior.WethenturntothecaseofBayesianLogisticRegressionunderthissameprior.We reviewthecannonicalapproachofperformingBayesianProbitRegressionthroughauxiliaryvariables, andextensionsofthistechniquetoBayesianLogisticRegressionandBayesianMultinomialRegression. Wethenturntothetaskoffeatureselection,outliningatrans-dimensionalMCMCapproachtovariable selectioninBayesianLogisticRegression.Finally,weturntothecaseofestimatingMAPparameters andperformingBayesianLogisticRegressionunderL1penaltiesandothersparsitypromotingpriors. 1Introduction Inthisproject,weexaminedthehighlypopularLogisticRegressionmodel.Thismodelhastradition- allybeenappealingduetoitsperformanceinclassication,thepotentialtouseitsoutputsasprobabilitic estimatessincetheyareintherange [0,1] ,andtheinterpretationofthecoefcientsintermsofthe’log- odds’ratio[1].Itisespeciallypopularinbiostatisticalapplicationswherebinaryclassicationtasks occurfrequently[1].Inthisrstpartofthereport,wereviewthismodel,itsmulti-classgeneralization, andstandardmethodsofperformingmaximumlikelihood(ML)ormaximumaposteriori(MAP)para- meterestimationunderazero-meanGaussianpriorfortheregressioncoefcients.Wethenturntothe caseofobtainingBayesianposteriordensityestimatesoftheregressioncoefcients.Inparticular,we examinetherecentlyproposedextensionsoftheBayesianProbitRegressionauxiliaryvariablemodelto theLogisticRegressionandMultinomialRegressionscenarios.Finally,weturntothechallengingtask

ofincorporatingfeatureselectionintothesemodels,focusingontrans-dimensionalsamplingmethods, andMAPand/orBayesianestimationunderpriorsthatencouragesparsity. 1.1BinaryLogisticRegressionModel Weuse X todenotethe n by p designmatrix,containing p featuresmeasuredfor n instances.Weuse y todenotethelength n classlabelvector,wherethevaluestakeoneither +1 or − 1 ,correspondingto theclasslabelforthe n th instance.Finally,wewilluse w torepresentthelength p vectorofparameters ofthemodel.Primarilyforeaseofpresentation,wewillnotaddressthe‘bias’term w 0 inthisdocument, butalltechniqueshereinareeasilymodiedtoincludeabiasterm.Underthestandard(binary)Logistic Regressionmodel,weexpresstheprobabilitythataninstance i belongstotheclass +1 as: 1 π(y i =+1 | x i ,w)= 1+exp( − w T x i ) (1) Forbinaryreponses,wecancomputetheprobabilityofthe’negative’classusingthesumruleof probability: π((y i = − 1 | x i ,w)=1 − π(y i =+1 | x i ,w) .WetypicallyassumeindependentGaussian priorswithmeansof0andvarianceof v onthecoefcientsofthemodel: w i ∼ N (0,v) )2( ToperformMAPparameterestimation,wetakethelogofthelikelihood(1)overallexamples,times theprior(2)overallparameters(ignoringthenormalizingconstant)togivethefollowingobjective function: n 1 X f= − log(1+exp( − y i w T x i )) − 2 w T w (3) v2 1=i Fromthisexpression,weseethattheMaximumLikelihoodestimateisobtainedbysetting v to ∞ . Differentiatingtheabovewithrespectto w toweobtainthefollowingexpressionsforthegradientand Hessian(using σ todenotethesigmoidfunction σ(x)=1/(1+exp( − x) ): 2

n w X g= − (1 − σ(y i w T x i ))y i x i − 2 v 1=i n 1 X H= − σ(w T x i )(1 − σ(w T x i ))x i x i T − 2 I p v 1=i ()4 )5( WenotethattheHessianisnegative-deniteandsubsequentlythattheoriginalfunctionis(log)concave, indicatingthatanylocalmaximizerofthisobjectivewillbeaglobalmaximizer.Asimplemethodto maximizethisobjectiveistorepeatNewtoniterationsstartingfromaninitialvalueof w untilthenorm ofthegradientissufcientlysmall(notingthatthatthegradientwillbe0atamaximizer).Thisresults inasimplexed-pointiterativeupdateasfollows: w=w+αH m − 1 g )6( Where H m isamodicationoftheHessiantobesufcientlynegative-denite,oranegative-denite approximationtothe(inverse)Hessian(see[2]).Thestepsize α canbesetto1,butconvergencemay behastenedbyusinglinesearchmethodssatisfyingsufcientdescentconditions(see[3]or[4]).We haveimplementedanapproachofthistypemakinguseofMatlab’s‘fminunc’functioninthedirectory ‘LOGREG’,andanexamplecallingthiscodeisincludedas‘exampleLOGREG.m’.Othermethodsfor optimizingthisobjectivearediscussedandcomparedin[4]. 1.2MultinomialLogisticRegressionModel ThebinaryLogisticRegressionmodelhasanaturalextensiontothecasewherethenumberofclasses, K ,isgreaterthan2.Thisisdoneusingthesoftmaxgeneralization[1]: exp(w k T x i ) π(y i,k | x i ,w)= P K exp(w T x i ) (7) j 1=j Inthiscase,wehaveamatrixoftargetlabels y thatis n by K ,and y(i,j) issetto +1 ifinstance i hasclass j ,and 0 otherwise.Theweightsarenowexpandedtoa p by K matrix,andwenowhave 3

anindividualweightvectorcorrespondingtoeachclass.Notethatwritingtheclassprobabilitiesinthis waymakesitclearthatweareemployinganexponentialfamilydistribution.Byobservingthatthe normalizingdenominatorenforcesthattheprobilitiessummedovertheclassesmustbeequalto 1 ,we cansettheparametervectorforoneoftheclassestobeavectorofzeros.Usingthis,wecanseethatin thiscasethesoftmaxlikelihoodwillbeidenticaltothebinarylogisticregressioncasewhenwehavetwo classes.Notealsothatthecoefcientsusedinasoftmaxfunctionretaintheirinterpretabilityintermsof changestothelog-odds,butthatthesechangesarenowrelativetotheclasswhoseparametersaresetto zero[1]. AgainassuminganindependentGaussianpriorontheelementsof w ,wecanwritethemulti-class penalizedlog-likelihoodforuseinMAPestimationasfollows: X n X K 1 X K f= − [y i w T x i − log(exp(w k T x i ))] − 2 w j T w j (8) i=1j=1 2v j=1 Abovewehaveintroduced y i asanindicatortoselecttheappropriatecolumnof w fortheinstace i . Weseethatthelog-likelihoodtermhasthefamiliar(numerator-denominator)form,subsequentlywe expectthegradientandHessiantocontainmomentsofthedistribution.IfweuseSM(i,k)todenotethe softmaxprobabilityofinstance i forclass k ,and δ i==j asthekroneckerdeltafunctionfor i and j ,we expressthegradientfortheparametersofclass k andtheHessianfortheparametersofclasses k and j asfollows: n 1 X g k = − [x i (y i − SM(i,k))] − 2 w k v 1=i n δ X H kj = − x i x i T [SM(i,k)(δ i=j − SM(i,j))] − j= 2 k v 1=i )9( )01( TheHessianremiainsnegativedenite,butnowhas (pK) 2 elementsinsteadof (pK) ,makingcomput- ingand/orinvertingtheHessianmuchmoreexpensive.Itisnoteworthythatinthesoftmaxcasewecan (andwill)haveahigherdegreeofcorrelationbetweenvariablesthanwedidforthebinarycasesincewe 4

haveadditionalcorrelationbetweentheclasses.WehaveimplementedMAPestimationformultinomial regressionmakinguseofMatlab’s‘fminunc’function(andhenceusingupdatesbasedonthegradient andaninverseHessianapproximationasdiscussedforthebinarycase)inthedirectory‘MLOGREG’, andanexamplecallingthiscodeisincludedas‘exampleMLOGREG.m’.Notethatthiscodeisnot vectorized,socouldbemademuchmoreefcient. 2BayesianAuxiliaryVariableMethods Abovewehavedescribedindetailthelogisticandmultinomialregressionmodels,andoverviewed somestraightforwardmethodstoperformMAPparameterestimationinsuchmodelsunderaGaussian prior.However,wewouldmuchratherbedoingBayesianparameterestimationinthesemodels,in ordertoobtainposteriordistributionsofthemodelparameters.WenowturntoBayesianmethods ofestimatingposteriordistributionsinlogisticregressionmodels.Inparticular,wewillfocusonthe Gibbssamplingmethodemployingauxiliaryvariablesandjointupdatestotheregressioncoefcients andauxiliaryvariablesproposedin[5]. 2.1BinaryProbitRegression Asdiscussedinclass,wecanderiveconjugatepriorsforthelogisticregressionlikelihoodfunction, buttheyarenotterriblyintuitive.Fortunately,wecantransformthemodelintoanequivalentformulation thatincludesauxiliaryvariables,andadmitsstandardconjugatepriorstothelikelihoodfunction.This methodisanextensionofthewell-knownauxiliaryvariablemethodforBinaryProbitRegressionof [6].Beforediscussingthelogisticregression,wewillrstreviewthesimplerBayesianBinaryProbit Regressionmodel,aspresentedin[5]. Using Φ todenotetheGaussiancumulativedistributionfunction,BinaryProbitRegressionusesthe followinglikelihood: π(y i =1 | x i )=Φ(x i T w) 5 )11(

SincethereisnoconjugatepriortotheGaussiancumulativedistributionfunction,weintroduceaset of n auxiliaryvariables z i with: z i =x i T w+ i )21( Here,  i ∼ N (0,1) ,and y i takesthevalueof 1 iff z i ispositive.Introducingtheauxiliaryvariables givesanequivalentmodel,butthismodelismoreamenabletosamplingsince w isremovedfromthe likelihood.IntheparticularcaseofaGaussianprioron w ,itadmitsastraightforwardGibbssampling strategywhere z i issampledfromindependent(univariate)truncatedGaussiandistributions,and w can canbesampledfromamultvariateGaussiandistribution.Specically,astraightforwardGibbssampling strategywith π(w) ∼ N (b,v) canbeimplementedusingthefollowing[5]: z i | w ∝ N (x i T w,1)I(z i >0)y i =1 T z i | w ∝ N (x i w,1)I(z i ≤ 0)y i 6 =0 w | z,y ∼ N (B,V) B=V(v − 1 b+X T z) V=(v − 1 +X T X) − 1 )31( )41( ()51 )61( )71( )81( Asbefore,wewillassumethatthemeanofthepriorontheregressioncoefcients b iszero.Follow- ingfromthis,wenoteabovethat B isthenormalequationsinLeastSquaresestimation,andthat V is thecorrespondinginverseHessian(correspondingtotheinverseconvarianceorprecisionmatrix).The Gibbssamplerproducedbytheabovestrategyistrivialtoimplement(givencurrentlyavailablelinear algebrasoftware),sinceitonlyrequiressamplingfromamultivariateGaussian,andtruncatedunivari- ateGaussians.Unfortunately,asdiscussedinclass,thisstraightforwardGibbssamplingapproachis innefcientsincetheelementof w arehighlycorrelatedwiththeelementsof z . 6

Tocombatthecorrelationinherentin w and z intheabovemodel,[5]proposedamethodtoupdate w and z jointlybyusingtheproductruletodecomposethejointprobabilityofthemodelasfollows: π(w,z | y)=π(z | y)π(w | z) )91( Theproposedmethodsampleseach z i fromaGaussiandistributionwithmeansandvariancesde- rivedfromaleave-one-outmarginalpredictivedensity(see(5)in[5]),updatingtheconditionalmeans aftereachupdateto z i ,thensampling w fromitsconditionalnormaldistributionafterallofthe z i have beensampled.Althoughresultsarepresentedshowingthatthisjointupdatingstrategyoffersanad- vantage,wewillnotreviewitindetailsincethereexistsasimplermethodtofacilitatejointupdating inlogisticregression.Nevertheless,wehaveimplementedthisblockupdatingscheme(basedonthe pseudocodepresentinthepaper)inthedirectory‘PROBREGSAMP’.Forourimplementations,we usedasimplerejectionsamplingapproachforsamplingfromtruncateddistributions,whereeachre- jectedsamplerestrictsthesamplingdensityenvelope.Anexampleshowinghowtorunthiscodeisat ‘examplePROBREGSAMP’. 2.2BinaryLogisticRegression BeginningfromthebinaryBayesianProbitRegressionmodelabove,[5]proposetoperformbinary BayesianLogisticRegressionbyreplacingtheindependentGaussianprioron  withindependentlogistic distributions.Unfortunately,thissignicantlycomplicatesthesimplesamplingstrategiesabove.To facilitatestraightforwardsamplingofthismodel[5]introduceanaddtionalsetofauxiliaryvariables λ 1:n andmodifythenoisefunctiontobeascalemixtureofnormalswithamarginallogisticdistribution asfollows(whereKSdenotestheKolmogorov-Smirnovdistribution):  i ∼ N (0,λ i ) λ i =(2ψ i ) 2 7 )02( 2()1

SK ∼ ψ i ()22 Ifwetemporarilyview λ asconstant,weseethatthisisidenticaltotheProbitmodelabove,except thateachvalueof z i hasanindividualterm λ i foritsnoisevariance.Subsequently,weknowhowto samplefromthismodelforxed λ .ItsimplyinvolvesusingWeightedLeastSquaresinsteadofLeast Squares(andassociatedinverseHessian),andsamplingfromindividualtruncatednormalsthathave differentvariances.SubsequentlywecanimplementaGibbssamplerifweareabletosamplefrom the KS distribution.Fortunately,[5]outlineanrejectionsamplingmethodtosimulatefromtheKS distributionusingtheGeneralizedInverseGaussianasthesamplingdensity.Thus,wecanimplement astraightforwardGibbssamplerforbinaryBayesianLogisticRegressionusingthisrejectionsampling approachinadditiontothefollowing: z i | w,λ ∝ N (x i T w,λ i )I(z i >0)ify i =1 z i | w,λ ∝ N (x i T w,λ i )I(z i ≤ 0)ify i 6 =0 w | z,y,λ ∼ N (B,V) B=V(v −1 b+X T Wz) V=(v − 1 +X T WX) − 1 W=diag(λ − 1 ) )32( )42( ()52 )62( )72( )82( )92( Twoapproachesarepresentedin[5]toperformblocksamplingoftheparameters.Therstuses thesamestrategyasintheProbitmodel,where w and z areupdatedjointlyinthesamewayusingthe abovemodicationstotheconditionals,followedbyanupdatetotheaddtionalauxiliaryvariables λ .We implementedthisstrategy(basedonthepseudocodefromthepaper)inthedirectory‘LOGREGSAMP’, ‘exampleLOGRESAMP’isanexamplescriptcallingthisfunction. 8

Thesecond(andtheauthor’spreferred)strategyforblocksamplingpresentedin[5]updates z and λ jointly,followedbyanupdateto w .Sampling w and λ remainsidenticalinthisapproach,butsampling z i becomeseasier.Inthisapproach, z i | w,λ followsatruncatedlogisticdistributionwithmean x i T w andascaleof 1 .Notonlydoesthisobviatetheneedforcomputingmarginalpredictivedensities, theinverseofthecumulativedistributionfunctionofthelogisticdistributionhasaclosedformand issubsequentlytrivialtosamplefrom(althoughweagainusedasimpleadaptiverejectionsampling techniqueinourimplementation).Weimplementedthisstrategy(basedonthepseudocodefromthe paper)inthedirectory‘LOGREGSAMP’,‘exampleLOGRESAMP2’isanexamplescriptcallingthis function. 2.3MultinomialLogisticRegression UnliketheProbitRegressioncase,thebinaryLogisticRegressionsamplingtechniquesabovehave atrivialextensiontothemulti-classscenario.Inadditiontohavinga y and w variableforeachclass aswesawinSection1,wenowhavea z and λ vectorforeachclass.TheGibbssamplerpresentedin [5]simplyloopsovertheclasses,performingthebinarylogisticregressionsamplingtechniqueforthe currentclasseskeepingallotherclassesxed.Weimplementedthisstrategy(basedonthepseudocode fromthepaper)inthedirectory‘MLOGREGSAMP’,‘exampleMLOGRESAMP’isanexamplescript callingthisfunction.Unfortunately,wefoundthatthisdoesnotmakeanespeciallyeffectivesampling strategy,andthatthetechniquestaysinareasofthedistributionthatwerefarfromtheMAPestimate, anddidnotproduceaccurateclassicationresults.Wehypothesizethatthisisduetoseveralfactors.The rstfactorissimplythelargernumberofparametersinthismodel.Anotherfactoristhat,asdiscussed previously,therecaninherentlybeamuchhigherdegreeofcorrelationinthesoftmaxcasethaninthe binaryscenario.Finally,wenotethatthesamplingstrategyofloopingovertheclasses,andrunningthe binarysamplerisnotespeciallycleveraboutdealingwiththesecorrelations,sinceitrequiresseperate samplingofthe z valuesforeachclassinadditiontothe w valuesforeachclass,andajointupdate wouldlikelyimprovetheperformance. Beforemovingontofeatureselection,Iwouldliketooutlinesomeextensionsoftheabovemodels 9

thatIwouldhavelikedtoexplored,ifIhadmoretime.Oneideawithsignicantpotentialforimprov- ingthesamplingstrategiesistointegrateoutparameters.Giventhehighdegreeofcorrelationbetween variables(especiallyinthemulti-classcase),thiswouldlikelyimprovethesamplingstrategiessignif- icantly.Anotherareaofexplorationistonottheviewcovarianceorhyper-parametersasxed,and exploreposteriorestimateswithpriorsonthesedistributions.Thisisespeciallyrelevantfromthepoint ofviewofmodelgeneralization,sincethecovarianceandhyper-parameterscansignicantlyaffectthe classicationperformanceofthemodel. 3FeatureSelection AmajorappealofLogisticRegression,besidesitsintuitivemulti-classgeneralization,istheinter- pretationofitscoefcients.Asdiscussedin[1],researchersoftenexploredifferentcombinationsofthe featuresinordertoproduceaparsimoniousregressionmodelthatstillprovideseffectivepredictionper- formance.Inthissection,wediscussautomatedapproachestothisfeatureselectionproblem.Werst presentanextensiontotheabovemodelsthatincorporatesfeatureselectionthroughtrans-dimensional sampling.Wethenturnourfocustopriorsthatencouragesparsityinthenalmodel. 3.1Trans-DimensionalSampling Focusingonthebinarylogisticregressionscenario,onemethodtoincorporatefeatureselectioninto theprocedureistoaddyetanothersetofauxiliaryvariables, γ 1:p .Specically,ifthebinaryvariable γ i issetto1thenthecorrespondingvariableisincludedinthemodel,andif γ i issetto0thenthe correspondingvariableisexcludedfromthemodel(ie.setto0).[5]proposesthismodel,andsug- gestusingthemodelpresentedearlierforbinarylogisticregressionwiththeseauxiliaryvariables,with jointupdatesto { z,λ } andto { γ,w } .Theyproposethat γ | z canbesampledusing(Reversible-Jump) Metropolis-Hastingssteps.Specically,sampling z , λ ,and w remainsthesame(butusingonlythe activecovariateset),andweacceptatrans-dimensionalstepfrom γ to γ ? (underasymmetricproposal) usingthefollowingacceptanceprobability: 01

| V γ ? | 1/2 | v γ | 1/2 exp(0.5B γ T V γ − 1 B γ ? ) α=min { 1, ??} (30) | V γ |1/2| v γ ? |1/2 exp(0.5B γ T V γ − 1 B γ ) [5]usesasimpleproposaldistribution,theyipthevalueofarandomlychosenelementof γ .We implementedsamplingfromtheabovemodelinthecaseofbinarylogisticregression(withfeatureselec- tion)inthedirectory‘LOGREGSAMPFS’,anexamplerunningthisroutineis‘exampleLOGREGSAMPFS’. 3.2PriorsEncouragingSparsity Althoughtheabovestrategytoincorporatefeatureselectionintothemodelisasimpleextensionofthe logisticregressionmodel,ithasmajordrawbacks.Specically,updatingsinglecomponentscausesvery slowexplorationofthespaceof 2 p variables.Asdiscussedinclass,wecouldjointlyupdatecorrelated componentstosignicantlyimprovetheresults.Analternatestrategy,especiallyrelevantwhen p isvery large,istousepriorsthatencouragesparsity. 3.3MAPEstimationoftheLogisticLASSO TheLASSOprioradvocatedin[7](bututlilizedearlierunderthename‘BasisFunctionPursuit’[8]) iscurrentlyapopularstrategyforenforcingsparsityintheweightsoftheregressioncoefcients.From thepointofviewofoptimization,theLASSOpriorconsistsofusingascaledvalueoftheL1-norm oftheweightsasthepenalty/regularizationfunction,insteadofthesquaredL2-normdiscussedearlier. Specically,ourobjectivefunctionbecomes: n X f= − log(1+exp( − y i w T x i )) − 1 || w || 1 (31) v 1=i Althoughtheaboveobjectiveisstillconcave,amajordisadvantageofthisobjectivefunctionisthat itisnon-differentiableatpointswhereany w i iszero.Hence,weneedtouseslightlylessgenericop- timizationapproachesforndingtheMAPestimates.Furthermore,wecannotuseefcientmethods suchastheonepresentedin[9]forLeastSquaresestimationunderanL1penaltyinordertooptimize thelogisticregressionlikelihoodfunction.Themostwidelyusedmethodforoptimizingthelogistic 11