DPMdencens之学解.docx
《DPMdencens之学解.docx》由会员分享，可在线阅读，更多相关《DPMdencens之学解.docx（24页珍藏版）》请在冰豆网上搜索。
DPMdencens之学解
DPMdencens{DPpackage}
RDocumentation
Bayesiandensityestimationforinterval-censoreddatausingaDPMofnormals
Description
ThisfunctiongeneratesaposteriordensitysampleforaDirichletprocessmixtureofnormalsmodelforinterval-censoreddata区间删失数据
.
Usage
DPMdencens（left,right,ngrid=100,grid=NULL,prior,mcmc,state,status）
Arguments
left
avectorormatrixgivingthelowerlimitforeachresponsevariable.NotethattheresponsesaredefinedontheentirereallineandthatunknownlimitsshouldbeindicatedbyNA.
right
avectorormatrixgivingtheupperlimitforeachresponsevariable.NotethattheresponsesaredefinedontheentirereallineandthatunknownlimitsshouldbeindicatedbyNA.
ngrid
numberofgridpointswherethedensityestimateisevaluated.Thedefaultvalueis100.
grid
matrixofdimensionngrid*nvarofgridpointswherethedensityestimateisevaluated.ThedefaultvalueisNULLandthegridischosenaccordingtotherangeoftheintervallimits.
prior
alistgivingthepriorinformation.Thelistincludesthefollowingparameter:
a0 and b0 givingthehyperparametersforpriordistributionoftheprecisionparameteroftheDirichletprocessprior, alpha givingthevalueoftheprecisionparameter（itmustbespecifiedif a0 ismissing,seedetailsbelow）, nu2 andpsiinv2 givingthehyperparametersoftheinvertedWishartpriordistributionforthescalematrix, Psi1,oftheinvertedWishartpartofthebaselinedistribution, tau1 and tau2 givingthehyperparametersforthegammapriordistributionofthescaleparameter k0 ofthenormalpartofthebaselinedistribution,m2 and s2 givingthemeanandthecovarianceofthenormalpriorforthemean, m1,ofthenormalcomponentofthebaselinedistribution,respectively, nu1 andpsiinv1 （itmustbespecifiedif nu2 ismissing,seedetailsbelow）givingthehyperparametersoftheinvertedWishartpartofthebaselinedistributionand, m1givingthemeanofthenormalpartofthebaselinedistribution（itmustbespecifiedif m2 ismissing,seedetailsbelow）and, k0 givingthescaleparameterofthenormalpartofthebaselinedistribution（itmustbespecifiedif tau1 ismissing,seedetailsbelow）.
mcmc
alistgivingtheMCMCparameters.Thelistmustincludethefollowingintegers:
nburn givingthenumberofburn-inscans, nskip givingthethinninginterval,nsave givingthetotalnumberofscanstobesaved,and ndisplay givingthenumberofsavedscanstobedisplayedonscreen（thefunctionreportsonthescreenwhenevery ndisplay iterationshavebeencarriedout）.
state
alistgivingthecurrentvalueoftheparameters.Thislistisusedifthecurrentanalysisisthecontinuationofapreviousanalysis.
status
alogicalvariableindicatingwhetherthisrunisnew（TRUE）orthecontinuationofapreviousanalysis（FALSE）.Inthelattercasethecurrentvalueoftheparametersmustbespecifiedintheobject state.
Details
ThisgenericfunctionfitsaDirichletprocessmixtureofnormalmodelfordensityestimation（EscobarandWest,1995）basedoninterval-censoreddata:
yijin[lij,uij）,i=1,…,n,j=1,…,m,
yi|mui,Sigmai~N（mui,Sigmai）,i=1,…,n,
（mui,Sigmai）|G~G,
G|alpha,G0~DP（alphaG0）,
where, yi=（yi1,…,yim）,andthebaselinedistributionistheconjugatenormal-inverted-Wishartdistribution,
G0=N（mu|m1,（1/k0）Sigma）IW（Sigma|nu1,psi1）
Tocompletethemodelspecification,independenthyperpriorsareassumed（optional）,
alpha|a0,b0~Gamma（a0,b0）
m1|m2,s2~N（m2,s2）
k0|tau1,tau2~Gamma（tau1/2,tau2/2）
psi1|nu2,psi2~IW（nu2,psi2）
Notethattheinverted-Wishartpriorisparametrizedsuchthatif A~IWq（nu,psi） then E（A）=psiinv/（nu-q-1）.
Toletpartofthebaselinedistributionfixedataparticularvalue,setthecorrespondinghyperparametersofthepriordistributionstoNULLinthehyperpriorspecificationofthemodel.
Althoughthebaselinedistribution, G0,isaconjugatepriorinthismodelspecification,analgorithmbasedonauxiliaryparametersisadopted.Specifically,thealgorithm8with m=1 ofNeal（2000）isconsideredinthe DPMdencens function.
Finally,notethatthisfunctioncanbeusedtofittheDPMofnormalsmodelforordinaldataproposedbyKottas,MuellerandQuintana（2005）.Inthiscase,thearbitrarycut-offpointsmustbespecifiedin left and right.Samplesfromthepredictivedistributioncontainedinthe（lastcolumns）oftheobjectrandsave（pleaseseebelow）canbeusedtoobtainanestimateofthecellprobabilities.
Value
Anobjectofclass DPMdencens representingtheDPmixtureofnormalsmodelfit.Genericfunctionssuchas print, summary,and plot havemethodstoshowtheresultsofthefit.Theresultsincludethebaselineparameters, alpha,andthenumberofclusters.
Thefunction DPrandom canbeusedtoextracttheposteriormeanofthesubject-specificmeansandcovariancematrices.
TheMCMCsamplesoftheparametersandtheerrorsinthemodelarestoredintheobject thetasave and randsave,respectively.Bothobjectsareincludedinthelistsave.state andarematriceswhichcanbeanalyzeddirectlybyfunctionsprovidedbythecodapackage.
Thelist state intheoutputobjectcontainsthecurrentvalueoftheparametersnecessarytorestarttheanalysis.Ifyouwanttospecifydifferentstartingvaluestorunmultiplechainsset status=TRUE andcreatetheliststatebasedonthisstartingvalues.Inthiscasethelist state mustincludethefollowingobjects:
ncluster
anintegergivingthenumberofclusters.
muclus
amatrixofdimension（nobservations+100）*（nvariables）givingthemeansoftheclusters（onlythefirst ncluster areconsideredtostartthechain）.
sigmaclus
amatrixofdimension（nobservations+100）*（（nvariables）*（（nvariables）+1）/2）givingthelowermatrixofthecovariancematrixoftheclusters（onlythefirstncluster areconsideredtostartthechain）.
ss
anintergervectordefiningtowhichofthe ncluster clusterseachobservationbelongs.
alpha
givingthevalueoftheprecisionparameter.
m1
givingthemeanofthenormalcomponentsofthebaselinedistribution.
k0
givingthescaleparameterofthenormalpartofthebaselinedistribution.
psi1
givingthescalematrixoftheinverted-Wishartpartofthebaselinedistribution.
y
givingthematrixofimputeddatapoints.
Author（s）
AlejandroJara
References
Escobar,M.D.andWest,M.（1995）BayesianDensityEstimationandInferenceUsingMixtures.JournaloftheAmericanStatisticalAssociation,90:
577-588.
Kottas,A.,Mueller,P.,Quintana,F.（2005）.NonparametricBayesianmodelingformultivariateordinaldata.JournalofComputationalandGraphicalStatistics,14:
610-625.
Neal,R.M.（2000）.MarkovChainsamplingmethodsforDirichletprocessmixturemodels.JournalofComputationalandGraphicalStatistics,9:
249-265.
SeeAlso
DPrandom, DPdensity
Examples
##Notrun:
####################################
#Bivariateexample:
#Censoreddataisartificially
#created
####################################
data（airquality）
缺失数据为NA
attach（airquality）
将数据中变量能直接读取
ozone<-Ozone**（1/3）
开立方根
radiation<-Solar.R
重新赋名
y<-na.omit（cbind（radiation,ozone））
删除带有na的行,并给出行号
#createcensored-data
xxlim<-seq（0,300,50）
yylim<-seq（1.5,5.5,1）
生成两个数列
left<-matrix（0,nrow=nrow（y）,ncol=2）
right<-matrix（0,nrow=nrow（y）,ncol=2）
生成与变量y同行,列为2的全0阵
for（iin1:
nrow（y））
{
left[i,1]<-NA
right[i,1]<-NA
if（y[i,1]for（jin1:
length（xxlim））
{
if（y[i,1]>=xxlim[j]）left[i,1]<-xxlim[j]
if（y[i,1]>=xxlim[j]）right[i,1]<-xxlim[j+1]
}
left[i,2]<-NA
right[i,2]<-NA
if（y[i,2]for（jin1:
length（yylim））
{
if（y[i,2]>=yylim[j]）left[i,2]<-yylim[j]
if（y[i,2]>=yylim[j]）right[i,2]<-yylim[j+1]
}
}
>left
[,1][,2]
[1,]1502.5
[2,]1002.5
[3,]1001.5
[4,]3002.5
[5,]2502.5
[6,]502.5
[7,]01.5
[8,]2502.5
[9,]2501.5
[10,]2501.5
[11,]502.5
[12,]3001.5
[13,]3002.5
[14,]501.5
[15,]3002.5
[16,]01.5
[17,]0NA
[18,]3001.5
[19,]01.5
[20,]502.5
[21,]02.5
[22,]2503.5
[23,]2004.5
[24,]2502.5
[25,]1002.5
[26,]2503.5
[27,]3002.5
[28,]1002.5
[29,]1502.5
[30,]2502.5
[31,]02.5
[32,]1001.5
[33,]1001.5
[34,]2504.5
[35,]2003.5
[36,]2002.5
[37,]1503.5
[38,]3002.5
[39,]2503.5
[40,]2504.5
[41,]2504.5
[42,]1503.5
[43,]2501.5
[44,]1502.5
[45,]01.5
[46,]2503.5
[47,]2502.5
[48,]2503.5
[49,]1503.5
[50,]2003.5
[51,]02.5
[52,]2503.5
[53,]2004.5
[54,]502.5
[55,]503.5
[56,]2003.5
[57,]2503.5
[58,]2503.5
[59,]2503.5
[60,]502.5
[61,]01.5
[62,]502.5
[63,]2504.5
[64,]2003.5
[65,]2004.5
[66,]1503.5
[67,]2502.5
[68,]1503.5
[69,]502.5
[70,]503.5
[71,]1002.5
[72,]2002.5
[73,]1503.5
[74,]2502.5
[75,]01.5
[76,]2003.5
[77,]2005.5
[78,]2003.5
[79,]2003.5
[80,]2004.5
[81,]2003.5
[82,]1503.5
[83,]1504.5
[84,]1503.5
[85,]1503.5
[86,]1503.5
[87,]503.5
[88,]502.5
[89,]2502.5
[90,]2002.5
[91,]2002.5
[92,]2502.5
[93,]2003.5
[94,]2502.5
[95,]2002.5
[96,]01.5
[97,]1001.5
[98,]2003.5
[99,]2002.5
[100,]01.5
[101,]2002.5
[102,]2002.5
[103,]2001.5
[104,]02.5
[105,]1002.5
[106,]01.5
[107,]01.5
[108,]1502.5
[109,]1501.5
[110,]1002.5
[111,]2002.5
>right
[,1][,2]
[1,]2003.5
[2,]1503.5
[3,]1502.5
[4,]NA3.5
[5,]3003.5
[6,]1003.5
[7,]502.5
[8,]3003.5
[9,]3002.5
[10,]3002.5
[11,]1003.5
[12,]NA2.5
[13,]NA3.5
升级会员 