第34卷第1期昆明理工大学学报{理工版}http://www. k田. 34 com/ 2009年2月Joumalof Kunming Univer百ityof Science and Technology (Science Feb. 2009 and Technology) 对带有有序分类变量的统计模型的贝叶斯分析付英姿，戴琳，叶凤英(昆明理工大学理学院，云南昆明白∞93) 摘要:针对带有有序分类变量的非线性再生散度因子分析模型建立起了一套贝叶斯方法，在计算过程中，采用了结合Gibbs抽样和M-H算法的混合MCMC算法来产生阀值，结构参数和潜在变量的联合贝叶斯估计，此外，还建立起了评估模型合理性的拟合优度统计量，贝叶斯方法的具体应用通过数据模拟加以说明.关键词:再生散度模型;Gibbs抽样;贝叶斯分析中固分类号: 文献标识码:A文章编号:1007-855X(2009)01 -0115 -06 Bayesian Analysis of Statistical Model with Continuous and Polytomous Variables FU Ying-zi, DAI Lin, YE Feng-ying (Faculty of Science, Kunming University of Science and Technology, Kunming 650093, China) Abstract: A factor analysis model with mixed continuous and polytomous variables is studied in this paper. A Bayesian approach is proposed to estimate simultaneously the thresholds, the structural parameters and potential variables via a hybrid Markov chain Monte Carlo algorithm that combines Gibbs sampler and Metropolis -Hasting algorithm. A goodness -of -fit statistics for assessing the plausibility of the posited model is introduced and the Bayesian procedure is illustrated by a simulation study. Key words: reproductive dispersion model; Gibbs sampler; Bayesian analysis o ~I言主要针对非线性再生散度因子分析模型建立起一套贝叶斯方法并且能够同时有效的处理混合RDM分布数据及有序分类数据.在此，采用数据添加的思想在后验分析中将观测数据集添加到模型中的潜在变量和有序分类数据中，同时在计算中采用Gibbs抽样技术川和MH算法[2]从后验分布中模拟出足够数量的观测并以此得到贝叶斯估计.除了点估计外，我们还建立起了获取标准差和拟合优度统计量的统计方法.1非线性再生散度因子分析模型根据Jorgensen[3]，如果随机变量Yi= (Yli’ ,Y.) T为px 1维的显变量，ι=(~(I)，…，~(q)) T是潜在p因子的随机向量，当￡给定时，假定yu，(i=1，1…n，k=I，2，…，p)有如下形式的概率密度函数:p(竹yL 2ψk J凡k阳"俨昨j 其中α(.，.):::0为某一适当的己知函数，d(.,. )为定义在Cxθ(θçCÇR，θ是一开区间，凸支撑点C为包含S的最小区间，S为概率密度函数的支撑点)上的单位偏差度函数(UnitDeviance)且满足条件:收稿日期:2∞8-05 -12.基金项目:昆明理工大学青年基金资助项目(项目编号:2007-82);省教育厅基金资助项目(项目编号:07WO∞1);国家社会科学基金资助项目(项目编号:07Bn∞1).第一作者简介:付英姿(1980-)，女，在读博士研究生.主要研究方向:应用统计.E-mail:fuyzI980@
116 昆明理工大学学报{理工版 (y; y) = 0, Vy eθ，d(y;μ) =O，Vy=i'μ其中μθ为位置参数或者表示随机变量y的均值，ψk> 0为散度参数，则我们称y属于再生散度模型(ReproductiveDispersion Model) ，简记为y-RDM(μ，ψ) .考虑以下因子分析模型:(KMR川队)μi = U + AF(ι) (1) k = 1, ,p, = 1, ,n, 其中μi= (μli'…，μpi) T ，u是px 1维的截距向量，A是pX r维的因子载荷矩阵，F(g)= (f1( ) , , J; ( ) ) T, (q < r) ’/1 ，.元，…d可以是这些潜在因子的非线性函数.类似于因子分析模型的一般假定，我们假定服从正态分布N(O，φ). 为了说明有序分类数据，不失一般性，令y=(yoT，yuT)T，其中YO=(Y1'仇，…，yJ是直接观测到的子集显变量且服从上述的RDM分布，而Yu= (y川，Yr+2，…，yp)是无法直接观测到的子集显变量，其信息、由观测得到的有序分类向量z= (z川，…，zp)T所提供，z与Yu的关系由一系列未知的阅值定义如下:当k= r + 1，…，p时r=C'α，"， < Yr+1 :5 a+,"d1 (2) zpαP句<Yp :5αP句+1其中Zk是属于0，1，…，b的整型变量，αkk= (叫，1，…，屿，b)T是阀值.在此，我们假定吨，0=∞，吨，b+1= kk+∞.对于=1，2，…爪，令Zι=(z川，川…，)T是相应于潜在变量Yui(Yr+1,i ，y吨H…，的观测，观测数据集为j(Y剧，z;)，= 1,2, ,nf. 2模型的贝叶斯分析令α=(α川，…，αp)，θ是包含了u，A，φ及ψ中所有未知参数的向量.令F= (˙1’ , n) T为潜在的因子矩阵，主要目的就是基于观测数据集|凡，zf建立起一套贝叶斯方法并对模型(1)作出相应的统计推断令p(θ，α)为θ和α的联合先验，由于z中的有序分类变量，似然函数p(凡，zIα，θ)涉及到一系列难以处理的多重积分，因此，直接处理后验密度p(凡，zIα，θ)是很困难的.根据数据添加的思想，我们视F与瓦为假定的缺失数据并在后验分析中把它们添加到观测数据集|凡，zf中进行分析，然而直接从p(θ，α|凡，瓦，F，z)中抽样仍然是一个困难，在此我们采用Gihhs抽样技术来从联合后验分布p(α，θ，F，瓦IY, oz)中抽取随机观测序列.我们以初值(α(0)θ(0),F(O) ，尺。))作为开始，以当前值(α(j)θ(J)，F(JUEiI))进行第j+ 1步迭代，抽样的程序如下:第1步，从p(FIθ(j)α(JUY;1)，凡，z)中抽取FU+I);第2步，从p(θF(j+1)α(JUY;J)，凡，z)中抽取θ(j+1 )第3步，从p(α，YuI F(J+1)θ(汁1)瓦，z)中抽取(α巾IUYY+I));重复以上步骤直至算法收敛-Geman[4l就曾经指出，在比较温和的条件下，对于足够大的J，(α(J)θ(J)，F(JUY;I)可以视为来自联合后验(α，θ，F，瓦|凡，z)的观测，于是后验分布就可以用相应的经验分布(αω，θ(J)，F(JUYY)U=J +1, l ,J + T)来近似，而算法的收敛性由相应于各个参数的EPSR(EstimatedPotential Scale Reduction[5)值所监控，如果所有未知参数的EPSR值都小于则算法收敛.条件分布在以下的部分中，我们简单的指出实施Gihhs抽样所需要的条件分布.考虑(θ，α，F，凡，瓦，z)的联合分布p(θ，α，F，Y，瓦，z)= p(θ，α)p(F Iθ，α)p(凡，凡，zI F，θ，α) o由F，Yo，瓦，z及α的定义易知:p(θ，α)= p(θ)p(α) ，p(FIθ，α) = p( F Iθ) ，且p(凡，瓦，zlθ，α，F)=
第1期忖英姿，戴琳，叶凤英:对带有有序分类变量的统计模型的贝叶斯分析117 p(YI F，θ)p(Z I瓦，α)，于是有:p(θ，α，F，Y'瓦，z)= p(θ)p(α)p(F Iθ)p(YI F，θ)p(z I瓦，α)(3) O我们首先考虑条件分布p(FIθ，α，孔，凡，z)，由(3)可得:p(FIθ，α，瓦，Yo，z)= TIp(ιI Yi，Zi'θ)民日(P(YiIι，Zi'θ)P(Çi Iθ) )于是:P(˙i I u,@,Y) oc exp{-过d(川kJ/l{lk+去(川k)-打了φ-IÇi}(4) 其中C(Yki'ψk)= log α(Yki，l{Ik)'我们接下来求Gibbs抽样第2步所需要的条件分布.根据Lindley和JSmith[6以及诸多先前对于NFA的贝叶斯分析的假定，我们考虑如下形式的共辄先验分布:p(u) -N(叫，立。)p(ψk -1) (αOAk ，ßo脑)Ak -N(A，ψkHOk)φ-IW(Ro,Po,q) ok其中JLT是的第k1于，!W表示的是逆Wishart分布，问，I, 0 ,UOAk , OAk ,AOk ,HOk ,HOk ,Ro ，Po分别是已被事先给定的超参数.以上形式的共辄先验分布在大多数的实际运用中具有很强的适应性，并且当数据容量大小合理时，超参数的选择很少影响到以后的分析.由(3)得:p(θI F, Y) oc p(θ)p(YI F，θ) 由于先验分布是相互独立的，在给定φ时，王的条件分布为N(O，φ)，同理可得:[φI FJ =IW(FFT+Ro-1,n+po,q) (5) 以ulMAKA)hxp{-÷z主d(川kJ/ψk才(u-u) T I, 0 -1 (u -u) } ( 6 ) oo州|川叫)叫k(l叫'Ak)exp{ -ßOA/ψk一过d(川kJ/仇+Zc(川k)} (7) 川AkI Y，u，ψ)仄叫-÷zd(川kJ/ψk-÷(AK-AOK)ν(A-A) } (8) kOk于是，Gibbs抽样第2步所需要的后验分布也得到了.最后考虑联合条件分布p(α，YI F，θ，凡，Z)，由于我们对阀值的选择信息所知甚少，为了处理更为普u遍的情形，我们采用以下的无信息先验:p(αk) = p(αυ，…，α-I) oc C αk，2 <…〈αk，bk-1，k=1,2, ,p 其中C为常数，同时，很自然的我们假定当k笋l时间，αt相互独立.对于k=r+l，…，p令Yuk = (Y, , klYkn) ,Zk = (Zkl’ ’kn) ，于是有:p(αk ，Y时，IF，θ，Y,z) = p(αk I勺，θ，F)p(Y叫，|αk，Zk'θ，F) O 其中p(αkI勺，θ，F)αTI1G ( Uk,Zki+) -G(问，Zμ)f ，且:1p( Yuk I问，勺，θ，F)仄TIexp~一土的μ)+ c(yψ) ~I(αk，Zki ,Uk,Zki+l) (Yki) (9) 问t-"r L 2ψk ~'Jki ’1-"kil ’ "\Jki ,’PkJ I G(, )为儿，的累积分布函数，!A(Y)为示性函数，相应的对于k=r+l，…，p，Y川儿川u时k川阳k们，YIZk'θ矶@川，盯盯F)仄H叫 丰2立};立d(y川付ρkiJ川+c( 川k忡冲ρ)}冲I净(αk札M川，阳吨Zk然后，和利l用多元MH抽样，α'y就能够从其边缘分布中进行有效的模拟.从而，Gibbs抽样第3步所需u要的条件分布也得到了.实施抽样从(5)中可以看出来，p(φIF)的分布是我们所熟悉的分布，从(5)中抽样是直接的并且是容易的.但是(4),(6) ,(7) ，(8)中的分布是非标准化的而且相当复杂，我们很难从直接其中抽取观测，在此，我们
118 昆明理工大学学报{理工版第34卷运用MH算法[2]解决这个困难.2为了执行MH算法，我们采用(4)中所给定的P(giIα，y，θ)作为目标函数，并选择N(O，σ0110)作为建议分布，41=φ-1+扫TA T I{I ~ ，其中.1=δaF(何ω5怡ιυJ饨|ι~i=O,I{I μpJIψ轧p)I ι~ιi严=0们，σ叫0J2为→参数，我们选择σ町02以使得5ιB平均接受率大致为.，则MH算法实施如下:以当2前值5ιω(j)进行第次迭代，然后从N川(0，σ町0110)中产生新的侯选值ι，接受侯选值5ιz的概率为:{lMi I Yi，Ø~J ’p(g,<j) I Yi'θ) 同理，我们继续采用MH算法来分别从密度函数p(ulY，F，Ak'制，p(ψkI Y,F,u,A) ,p(AI Y，u，ψ) kk 中抽取观测，其各自的建议分布分别为:N(O，σ12(1)，N(O，σ/lh)及N(O，σ32~)，其接受的概率分别如下:k{lP(u I Y,F ,Ak ，I{ILlminJt州kI Y,F,u,AL lminI |队ψ;，Jω'p(ttJ)|Y，F，AK，ψ)J L-'p(ψk(j) I Y,F,u,Ak)J L-'p(AI Y，u，ψ) J k最后，我们采用多变量的MH算法来从p(屿'YI勺，θ，F)，k= r+l，…，p中抽取测，第j步的MH算uk法的实施步骤如下:首先我们从以下的截尾正态中抽取侯选的阀值向量(αυ，…，屿，bk-l): 2αk，z -N(α己，σak) I(α-l ，αJILl)(α)， z = 2, ,b-1 (11) k 这里α已是Gibbs抽样在第j次迭代时间z的当前值，σJ的选择是为了获得接近的接受率，根据(αk(J+l)Cowls[7]的论证，我们接受侯选向量(αk'YYdW))此)作为新的观测i，的概率为minj1 Rk ，其中,/ 丁αk(J(j) αF)p(此i|勺，θukIk'Y，αYθF)'Yk>峙，，Zk'，--Ei-R, -= k (j)Yd(J)p(αk(JUYI时勺，θp(，F)αYIαkω，，勺，θF)k'uk ，根据(9)0)，(1可得:φ1(αJILI-d)/σαJ -φiWUI-dJWhiG(ak，Zki+1)-G(αâ) Rk =旨白(j) \ / _ J X 11 f’ / (j) f’ I (i) r (12 ) M2φiCαk，z+lαk，.)1σαk I -φiCαkμ-αk，.)1σαk lHG(αJf时1)-G(α(怡，估计及拟合优度评估令1(α(叫，θ(m)，F(叫，Y(m) ) : m = 1，2，…，M/是从恼，θ，F，YI }年，Z/中产生的随机观测序列，则α，uuθ，F的联合贝叶斯估计由其观测序列的样本均值得到:酒MM 1 = M-=M-IEFM) Tra=MIZFWOF 三卢显然，贝叶斯估计是其后验均值的一致估计，可是要推导出协差阵Vαr(α|Z)VIαr(θ|凡Z)，凡，Var(gi |凡Z)，的解析形式是十分困难的，基于模拟的观测，我们仍然可以通过样本协差阵去估计总体协差阵，例如:M _-l叫(a(m) α(Vα叫αIz) Yo(M- )T -â)(,= l)(13 ) I,此外，αθF，，标准差则可以通过以上矩阵的对角元得到.对假定模型的合理性评估一向是数据分析的最基本的问题.在贝叶斯的框架下，GelmanMeng和Sternl8J提出了一个后验预测性p值来作为评估模型合理性的拟合优度统计量ppp-，值被定义为:REP yD( IθFPB = Pr 1 ，，YJ三D(YIθF，，瓦)I YzHO ,,o I ~EP其中表示Y的一个复制D(.，1. )表示某个偏差变量.在此，我们选择如下的偏差变量:EP yRD( I θFY) (Y?EP-μì)，， 了。-川、飞Yu=工EP 其中ψdiag(R=ψ1 p) , 显然D(yIθF,l{I，，，，瓦)服从i(np)，于是ppp-值又等于:2(Yz) rPB( ,州叩)三D(YIθFY))FY，，州u，Io =，凡，川剧FdY户uu
第1期付英姿，戴琳，叶凤英:对带有有序分类变量的统计模型的贝叶新分析119 它的一个Rao-Bladkwelized类型的估计就为:1 2PB(YO’Z) = r-Lprob(x(np)三D(Yo，YJt)|θ(t),F(t) )) (14 ) 显然，PB(凡，Z)的计算是直接的.如果PB(凡，z，r)和的差距不大，比如在和之间，那么，我们认为(1)中假定的模型是合理的.3数据模拟运用数据模拟的结果来研究贝叶斯方法的经验表现采用logGamma分布作为RDM分布的特例.首先，由(1)中所定义的服从logGamma分布的因子分析模型中产生观测数据集|YKJ=1，…爪k= 1,2, 1 ,7 f ，即Yki-exp lYki -μki -exp( Yki -μkJ f ，其中，α(Yki，I/!k)= e-，ψk = ，d(Yki'的J= -21Yki-μki -exp(Yki -ι) + 1 f ，并且向=U+ A/F(~J ，A/表示A的第k行.此外，模型(1)中由潜在变量(ι1, k ~i2 )所构成的非线性函数F(ι)= (~i1'~泣'~i1~i2'~i2~i2) T ，连续观测Y7i又通过阅值(-1.γ， ,, 1. O’ )转换成有序分类观测气，其中-1. 2'和'被视为确定的已知值，于是，最终的观测数据集就为(Y1i ’Y2i’ ’Y6i ,Zi) ，其中前6个变量为服从logGamma分布的连续型变量，而最后一个变量为有序分类变(1λ21 0 0 0 0 0飞T10 0 1λ'" 0 0 0 1 量.载荷矩阵A为:A= 1- -~ , 4L ---1 10 0 0 0 1λ63 0 1 \o 0 0 0 0 0 λ74 ) A中的0和1被视为确定的已知值，Àij为待估参数.模型(1)中所有的未知参数的真实值分别为λ21-λ42 =λ63 =λ74 = , (φ11 ，φ12 ，φ22) = (1. 0, O. 5 ,1. 0) ，叫=O. 36 , (i = 1，…，7).在此因子分析模型中，一共有16个未知参数(U1，…，屿，λ21，λ42 ，λ63 ，λ74，φ11'φ12 ，φ22，αk1 ，αk2) ，样本量n为500.采用3组不同的初值来计算以上未知参数的贝叶斯估计，同时，先验分布中的超参数的值给定如下:1U o ,AOk f等于其真实值，立。=鸟，αOAk= ，βω= ,Po = ,Ro = 5矶，HOk为对角元都为的对角矩阵，以上这种情况可以被视为具有良好先验信息的情形.运用()所介绍的25 Gibbs抽样和MH算法来20 获得未知参数的贝叶斯估计.此外，为了降低随机误Iil 15 差的影响，我们将整个计5于算过程重复100次.在建议10 分布中，分别令σ02=2 ，σ2 = ，σ3= 5 2 和σak= 以使得平均1. 0接受率分别为，，137 4ο9 681 95312251497176920412313 2585 2857 3129 34013673 3945 和值见图1模拟研究中的EPSR值图1，发现当迭代次数大于 The EPSR value in the simulation study 1500次以后，所有参数的EPSR值于，保守起见，我们收集了2∞0次以后的总数为M= 3∞0个观测，根据(19)就能够得到16个未知参数的贝叶斯估计.估计的均值，估计的偏差(即未知参数的真实值与100次复制所产生的贝叶斯估计的均值之间的差异)及估计的标准差(SD)在表(1)中有所体现.此外，通过计算，模型的ppp-值等于，这说明模型和数据拟合较好.
120 昆明理工大学学报(理工版) 第34卷表1数据模拟中所有未知参数的贝时斯估计Tab. 1 The Bayesian estimates and the standard errors in the simulation study Parα EST Biω SD Para EST Bias SD UO. 368 λ21 1 λ42 2 6 λ63 O. 615 O. 015 3 O. 023 4 -0 4 O.∞70 UO. 355 ∞5 6 φ1I s UO. 372 φ12 O. 555 0,0599 6 U, O. 347 2 αkl 0,017 9 αn 参考文献:[ 1 J Metropolis N, Rosenhluth A W, Rosenhluth M N, et al. Equations of State Calculations hy Fastcomputing Machine [J J, Jour›nal of Chemical Physics. 1953,21,1087 -1091. [2 J HASTINGS W K. Monte C町10Sampling Methodsusi吨MarkovChains and Their Application [ J J. Biomet出,57 :97 -109. [3 J JORGENSEN B. The Theoηof Dispersion Models [ M J. Chapman and Hall, London, 1997, [4 J GELMAN S, GEMAN D. Stochastic Relaxation, Gihhs Di刨出ution，and the Bayesian Restoration of Images [ J J. IEEE Trans›actions on Pattern Analysis and Machine Intelligence3, 1984 ,6: 721 -741. [5J GELMAN A. Markov chain Monte Carlo in Practice[MJ. London: Chapman and Hall,1996. [6 J LINDLEY D V, SMITH A F M. Bayes Estimates for the Linear Model (with Discussion[ JJ. Journal of the Royal Statistical So›ciety, 1972, Series B ,34: 1 -42. [7 J COWLES M K, CARLIN B P. Markov Chain Monte Carlo Convergence Diagnostics[ J]. A Comparative Review. Journal of the American Statistical Association, 1996,91,883 -904. [8 J GELMAN A, MENG X L, STERN H. Posterior Predictive Assessment of Model Fitness via Realized Discrepancie恐[JJ. Statis›tical Sinica, 1996,6:733 -807. [9J吴刘仓，李会琼.兰辛森林数据的空间统计分析[JJ，昆明理工大学学报:理工版，2008，33(2):112 -117. 唱加甘情甘甜也乒昏世串回串哩庐喃晶咂府甘骨哩""哩府叩晶q晶q骨喃串哩府也/昏暗局审毒也严岛也乒昏咀庐飞g串也卢哩局也".，.也卢哩局也严Õ"，哩卢也乒霍、~~周也卢响周恩、可a>，、电府、唱属》串、啕品、咽，串叩'""，由南-局(上接第110页)3结论当然我们还可以从数学的立场对前面的公式进行更为详细的分析.是因为我们的目的是借助"数学模型"这个工具，来探讨"cp-城市吸引力如果一个城市没有吸引力，谈何发展.城市之间的竞争，很大程度就是人才的竞争[4-5] 最后，在我们的模型中，假定流进和流出的每个人都有相同的统计权重.若以对社会的贡献而言，不争的事实是有一般人才和特殊人才的区别.由于原始数据的缺失，无法补充例证.但即使是这样，也不失为一种有益的探索和尝试.当然也期待下一步能有实际例证来验证这个设想，同时也期望能从经济、金融、吞吐量等等多种角度来解说"cp-城市吸引力参考文献:[ 1 J王洋李翠霞.西方人口流动理论经典模型分析[JJ.东北农业大学学报:社会科学版，2006，4(3):74 -75. [2J刘渝琳，李洁.西部大开发中人才"回流"的研究[JJ.林业调查规划，2∞1，26(3):25 -29. [3J张海鹏，西方人口流动模型及其对我国就业政策的启示[JJ .山西高等学校社会科学学报，2∞5，16(1):76 -78. [4J王桂芝，袁博.有序人口流动模型及其实证分析[JJ.安徽农业科学，2∞7，35(29) :9367 -9369. [5J赖小琼.中国转型时期的人口流动[JJ.中国经济问题，2∞7，23(1):41 -47.
