﻿
 UFDC Home myUFDC Home  |   Help
<%BANNER%>

# Theoretical and Methodological Developments for Markov Chain Monte Carlo Algorithms for Bayesian Regression

## Material Information

Title: Theoretical and Methodological Developments for Markov Chain Monte Carlo Algorithms for Bayesian Regression
Physical Description: 1 online resource (94 p.)
Language: english
Creator: Roy, Vivekananda
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2008

## Subjects

Subjects / Keywords: bayesian, da, efficiency, markov, monte, multivariate, probit, px, regenerative
Statistics -- Dissertations, Academic -- UF
Genre: Statistics thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

## Notes

Abstract: I develop theoretical and methodological results for Markov chain Monte Carlo (MCMC) algorithms for two different Bayesian regression models. First, I consider a probit regression problem in which $Y_1,\dots,Y_n$ are independent Bernoulli random variables such that $\Pr(Y_i =1) = \Phi(x_i^T \beta)$ where $x_i$ is a $p$-dimensional vector of known covariates associated with $Y_i$, $\beta$ is a $p$-dimensional vector of unknown regression coefficients and $\Phi(\cdot)$ denotes the standard normal distribution function. I study two frequently used MCMC algorithms for exploring the intractable posterior density that results when the probit regression likelihood is combined with a flat prior on $\beta$. These algorithms are Albert and Chib's data augmentation algorithm and Liu and Wu's PX-DA algorithm. I prove that both of these algorithms converge at a geometric rate, which ensures the existence of central limit theorems (CLTs) for ergodic averages under a second moment condition. While these two algorithms are essentially equivalent in terms of computational complexity, I show that the PX-DA algorithm is theoretically more efficient in the sense that the asymptotic variance in the CLT under the PX-DA algorithm is no larger than that under Albert and Chib's algorithm. A simple, consistent estimator of the asymptotic variance in the CLT is constructed using regeneration. As an illustration, I apply my results to van Dyk and Meng's lupus data. In this particular example, the estimated asymptotic relative efficiency of the PX-DA algorithm with respect to Albert and Chib's algorithm is about 65, which demonstrates that huge gains in efficiency are possible by using PX-DA. Second, I consider multivariate regression models where the distribution of the errors is a scale mixture of normals. Let $\pi$ denote the posterior density that results when the likelihood of $n$ observations from the corresponding regression model is combined with the standard non-informative prior. I provide necessary and sufficient condition for the propriety of the posterior distribution, $\pi$. I develop two MCMC algorithms that can be used to explore the intractable density $\pi$. These algorithms are the data augmentation algorithm and the Haar PX-DA algorithm. I compare the two algorithms in terms of efficiency ordering. I establish drift and minorization conditions to study the convergence rates of these algorithms.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Vivekananda Roy.
Thesis: Thesis (Ph.D.)--University of Florida, 2008.
Electronic Access: RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2010-08-31

## Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2008
System ID: UFE0022377:00001

## Material Information

Title: Theoretical and Methodological Developments for Markov Chain Monte Carlo Algorithms for Bayesian Regression
Physical Description: 1 online resource (94 p.)
Language: english
Creator: Roy, Vivekananda
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2008

## Subjects

Subjects / Keywords: bayesian, da, efficiency, markov, monte, multivariate, probit, px, regenerative
Statistics -- Dissertations, Academic -- UF
Genre: Statistics thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

## Notes

Abstract: I develop theoretical and methodological results for Markov chain Monte Carlo (MCMC) algorithms for two different Bayesian regression models. First, I consider a probit regression problem in which $Y_1,\dots,Y_n$ are independent Bernoulli random variables such that $\Pr(Y_i =1) = \Phi(x_i^T \beta)$ where $x_i$ is a $p$-dimensional vector of known covariates associated with $Y_i$, $\beta$ is a $p$-dimensional vector of unknown regression coefficients and $\Phi(\cdot)$ denotes the standard normal distribution function. I study two frequently used MCMC algorithms for exploring the intractable posterior density that results when the probit regression likelihood is combined with a flat prior on $\beta$. These algorithms are Albert and Chib's data augmentation algorithm and Liu and Wu's PX-DA algorithm. I prove that both of these algorithms converge at a geometric rate, which ensures the existence of central limit theorems (CLTs) for ergodic averages under a second moment condition. While these two algorithms are essentially equivalent in terms of computational complexity, I show that the PX-DA algorithm is theoretically more efficient in the sense that the asymptotic variance in the CLT under the PX-DA algorithm is no larger than that under Albert and Chib's algorithm. A simple, consistent estimator of the asymptotic variance in the CLT is constructed using regeneration. As an illustration, I apply my results to van Dyk and Meng's lupus data. In this particular example, the estimated asymptotic relative efficiency of the PX-DA algorithm with respect to Albert and Chib's algorithm is about 65, which demonstrates that huge gains in efficiency are possible by using PX-DA. Second, I consider multivariate regression models where the distribution of the errors is a scale mixture of normals. Let $\pi$ denote the posterior density that results when the likelihood of $n$ observations from the corresponding regression model is combined with the standard non-informative prior. I provide necessary and sufficient condition for the propriety of the posterior distribution, $\pi$. I develop two MCMC algorithms that can be used to explore the intractable density $\pi$. These algorithms are the data augmentation algorithm and the Haar PX-DA algorithm. I compare the two algorithms in terms of efficiency ordering. I establish drift and minorization conditions to study the convergence rates of these algorithms.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Vivekananda Roy.
Thesis: Thesis (Ph.D.)--University of Florida, 2008.
Electronic Access: RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2010-08-31

## Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2008
System ID: UFE0022377:00001

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E20101221_AAAAAE INGEST_TIME 2010-12-21T06:50:03Z PACKAGE UFE0022377_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
FILE SIZE 32630 DFID F20101221_AAADHU ORIGIN DEPOSITOR PATH roy_v_Page_09.QC.jpg GLOBAL false PRESERVATION BIT MESSAGE_DIGEST ALGORITHM MD5
a1212bb2daf0ef429e799cc3a97c00c7
SHA-1
1052f9960fe87b0ee480ccef2d983be5d42afb9c
6eebc8ccfd9d12a74a20176fa823499a
4feca1ed2d6f72e3a314f0a2cb33cd771b7573f6
9c27d62f1f0c9865f6ac72d458e3ee15cf9c3c5f
b93092f0865b09c11a59a7a4a8b5ba17787b7f65
8efa872404e6f7cf1bc31aa098b552e2
4044bd00c7798c66c2037887f949589bc3b12f01
f758892c9bc9a39c92aeddbae667586b
37c024a987823fdd08e1cb24d759f028ce348320
2a85af5fa66222295c0cf60ab1affd34
2510b931ea564309136c2216905cf09c4384ed89
75af7de1e93abc6e4350864f5fca240a
1b5707db7323274690c9e65258137737d3eb2e7b
554a5a8ecbf92f45e180f1be1c888b46
ac4f1260480ac22be4877e400e7aa530ac4d9b4f
3a7ac8da5ab15fc04eca999b1f636b51
226a78a8ba256778fe3a616f824f9ac917f9e611
5df53ae7ec5f3c5ff3e6960248b4220c
70920abf4f72bd40f2753ee6db26ac805ac9e792
8f78226835b28bf78f2765b8016a59e63ba90e24
b97049522d3ec6fc92924924d65f35a1
bf619885f2b5d721e2cabfdb7b7e6107b92c7c5b
a73313c32bfb30195f97f1e71a35b7d5
63592e3803c5ec26b7fc2c5df729258d
1b4ff5cd483650d9666c7c81d2fca129659f8983
d7dce1dc2a83faff6b399afed7612a25
b743e5cb2c2ce7b66521657ab3a77595ccf53550
5b9c9db7668d9d30daced6261abe8e66165f697e
5487fc825e5f96d03d8b934d519f8130e2fb0c11
ec220384b76e944ff97f3e5b9d10285f
38f737fbb9f28bc0c65994317c387b57be08a3f2
2c34674ff819ae775acbecd50bdf3ea3
398ba1eea2f9b4a46ecc1350ebf5c2b6b425d0b3
ee7df206dcba369624487e7564838eb7
8d80f392b52c5c7d02b90f96a37699e7de6b81ee
f32bfbb1ce31eebab2620e91239df365
9dd99ae1d6559db6567ecfb1cfe4d88f3d2bcc55
cfd97f2b9666fd1c3a5a4b9c9843ca5f
127ee16218c71c5dd253da56e43c0b7af1e37808
ce1b065bcc1ff5d80f1dd13120e7245627c6bb65
12c4a798a94c3101613a69b9d0212c96
fb951fd4b678b47aa852674e4256e4f0e9b50ec4
c8f0da36d4348807a0e165c0d7e9f007
7bccec8458f573795045b2fe22e5777d3999233b
7bb57f8c39cb693c63dec3f1e0ddff20
7e175cd9be8996f4b50b9fa50823aea8dde6df80
8f8cf0d65acc99529225a39481f227ba
252b6d4be600c67559323b5deefca1595bbf9004
36965c13f42656dd753b1e945f77332d
1706c5e3a991f083970f63ee9b4df8c6e35efc60
6384ba1f1b64d2673a862ded8e2aaf48
56257d710b1f068518620a918426bf8b75a1b807
f72ff715c90f3198367888582950c602
cdd2d3e021c598165ff88502d2415c4030dcddc8
ac30e080f00c9fdc5f82e6c4f56a73e3
f123d95695acdb93982cdc2f9ee7c2a5a5c6614b
c63a1c50610d60ec70e72815103592eb
d88c9668a23a47723327151e858c81656e9364e7
8e92de840be9d9ebde3f389eaed1f138
c37fa23aa6639e89d868d18ee24431ac18c68dae
22d06bdc8d67b4f732186afbc0b00f87
1444728af20771cf9267b6173a7c6bddf2e13f1e
afe324191e790ac28f2050094c47cba7
01d68e35876a8d215bf6e643f82a42d99188ecc9
4c105ce90baaebde379ca02ac15e13b3
5aa5f0381935fce39cd835b19a0cb4603db80c16
b83210a4a592564298ea4d92fc1f518f
e0027ae1a753ee22228cf0f19d7b2a831d7f3506
90e60062b12773929f573c2338a98994
e6f8a005ed547695b070bd38e9207a49a33b74c5
6182ef38423e5e8bed9026eff1ea1117
372bb8ee271f1f48706f538174ec3108
11ed256b38d49ea4a109a7e7bea8ace82dde0fdc
424e20a412a43f8bc16bd10fca3befe6
8ff35a242d2b596b5f80d673b54b5f2544b7e139
4f1e1d4a7577309290f57784e057252b
fb3fe3775137b66d8f6163ec35a51bc255ab61ee
d298fa94fb8ddf8f782af738daea1d4c346156b1
9f47f119f039853113523912619f262f
63ce525e19388dbe1eea68b11631baf424610f0d
39026c98e49a502f4522a986ec4cb19c
0ec4f57d76b521abc86ebb9b3a41bd551ca69463
0dff7cd788764632fe3f4d9bdb6399cb27c36fa4
72db2e240bf83016f424be53d1c4bee8
b88499a0631212e39f922224b04cf9f9e3e2d494
aaa338141198a2541f08d4b3f9b65d80
c640208c1b727330674925bd0fb094450069ac77
c65afabf4877eb49b761a889621ca723
99c4752ee5776ff24f47924209721f8a97288b80
4d90f751035db282ce40a4534414363a
7cc4ef1ede1c5f19689d7ec7774d02cde953fa41
b5bbd4facb86d0f8243b45311b2cb5ac
1bc1c06388db4674f653d44471b3960e6f8e8fc4
084887317b7cf0eda175bc943b8c574c
d2e6d6709d683041999e140f8623cf57620e20c0
f607be0f851d547b139c28a5805d1942
8093478bba91d279f28f7b2fe8b9c7a4e7631dac
286d1471cd586fea6d9c0a0a1800c2b3
7606a343be54fbc7e5a8c216e63092b2
c494364d154e1898cec2b2c4f388017b6e7f2482
c78de6fdcbaa6189155b8dce87bb4d52
b3cf68330134673a2bccef4a524f51fce923cea9
72ce87bd27fe3c0100b663080cedbd44
6b22b2ee48317be867b75c00533f1c5ed1fdd55d
e74a462546a9d6297e0c8818a42ee44b
301335f988864621b1913fb142cc7c53b4a1ccbb
4750f8d664a1e198ceb78da21f4a20d6
5858778672ac813069575650fec3fa678919cb46
1d5c7b908138ab380ba9505f145151b3
938de4b8843ac982d33a1bef4f6a244a7c1a8e6f
a67729f089332af51a161e5e68d51c31
df0aef9cff1cb56133977aa355eca54aa9b26ca7
b8dea69d2938441c3eb07d9ae972065b
6b19e1b152f53ba1523288d60e6f68a0324388ec
f745a713b4df682b53881d6c3c4ecf6c
58afb85759d47da8a38b83bee1e13f00
2d3cc4c690ee79ec27af903e8612c7607ea8cd21
443dbfa8980717708ca785e2d9fab692
b7e347c3158331dac4d4d9791f1f54af1f48f47d
9be8a5b1e2bd736e14f2403fc8b91f05
296928150b5b07d57b3f75522b267318
6bccf2efdc308c6311dbac0e3c09a38c8f56fbb6
f326bdb746c2538ca3f35506e12f9809
81de06d7e804b4f1309a394ae45a2bee445b87fb
0843391f8a4239ac777648c65f4e4f2c
96ee77964e81fed289307b3307bf70f301ed0bdf
77411dbd93e07e1ac5afe75c13f4c911
7307d2e4aab16657e7d0c00bfa931cf0cf876a66
68e805cbd0dcfabec3c3b2227fd6bab7
7cf4f8a22cddf084d3eb6f732eeb73e73d4dd928
ecfd7d9df463e5f06c2fe9c53c9ab631
4136919a74b09e52a22c7a733f0711663283719e
10b1cb81e9c4a7d80308ec7aa6be38a8
6eac96d0fe298330e7826f255e1e4e54726c3c5d
0ec4a9fc6d2c3a22aec4399aafdfdacd
3351cf11ae0e180e1a8d4ea9d54c6f11dd2f9731
8a3061eaebd3a52af9678a01395c763c
8277a79bfb75f57bd93b1bf52399429190ef2b0f
24dce7317e6b8925893778c3a279ffd3
2e4abaa7700f9dfef57b8a1afe48b29034acbccf
da2e1358e94c21ed92defd8d97e276a8
53d7bd708140b0c6defe330e9908b4b02354fe03
41ffce568127754b599722951b6e34a5
4e7d1fa81bd6f7c9d7c464f8233172bdacfc66d4
a8ee286533bc3c0bd585a046e5746350
2e41e4df7794a6dbfb1e1b6820e32f123eff8b35
bcd77841b91845b936a70a26b60b62c7
2b6f39d570d541d232412b3a374d297176b75d60
33fc8c5322b43d26edf2efcb33069189
724d2a40ddf440b8d1af9625b29101abb5d84798
cf70c1e6a2b55a34a793de1526b4e837
fbe75d2a3bfd673ded739cd9a37e428344318ab3
3f44a8feeaf5bebdb5421e767c7c07b0
e1bf6a3d45a1d88531aff465be9290f6a631fc89
b5b7cc4dbba0cc0b171df97065aeabdf
3bb051351601c0ca96a4dd9159aed98bb7c97248
c2cb9f79a4741193c36a142b1f7f13d4
c8160ee785bef21f385a91190c77a37540337767
80c2e615a5eda343e2887e41531c677f
ed96bc158b52abba4bc7f22a05e0e34b2723ae22
81de12c4f62024c4508a7103f9cd9037
bacd62ec414d2894edfda23eab1968142c20d513
6805d9683eee0b8e6880e18bd3d3e38f
3e87feeb497ec8ec829dbc566979d33020c92714
3cf7491f31314d60ce72162bd80ff53a
1c9776e674beda277f3bd26b1f0af2324af3fae6
1f621e6dc42fde1bb59321cff936918f
8204c16ac172fd11387395ae140f64542a42c9c4
e40eef7ae3ef6055e45b16592e43a5f9
b79571848850dca03d5e3859ccede6e1fbf55f3e
fdec4d1805171310c12cf0de4ab38539
e9253930240e37160ed86824e01b7655
570f76260c7da8625d835d6793839cdb01b9876f
f7971a08ece2ee2e06b7eb9c90836e26
6fa116617ee3ec4b8d44350b6c607152a11e07dc
7a82b9ebe80b5dbbbb7dc8f767966625
349b7fb21651ac7dd270b9421f8d6ff6ab5a05d4
ddd55310c79df8772944d672570948f5
70f1767acbea1f92c2992bf8e431932df4f1db9b
375c24ed97187cb964362875b9f1f65d
fcd62ba32e4bf96d63b9e6fd72d7eb170fa4449b
f7ec90c42c0c13481a3d432f30f8de36
8209316cd514e43ac1700b9def6a4e0e
d45b9c9277ca6cf8a8c6a9b6959249e1d5126edb
351d8237dcb9688820922460560ecfe4
62c208600637dba00b406bfcab7587247e2ffb7a
8f4928858dd765d7f32bd40682685b72
ac3c49013944ab9c75fe1f052487f867bfed05a6
52d288f46765c0535479100a7f9850aef0d778b1
b82e63bd013e3cf44025b19019da628a
24de80d0f60be4b70c0c203c41f40fc8c28307b5
8ec64ca2b3176a22a6c1d55ace4d8033
208d1593d50940db3b72dacf302dd84aa58e5ab9
573f762be4601dfe3200f267726ed05c
f1db42cbbd6bbe2503fda012da7c6d5af3ec8251
bb23bb8959dcc4d9959131c45b5b0122
2a5b8e9605d17a3233ff8560e70ae164
2b3ea6e305ac37281b47d4aec7c1bc26854f213b
ced69128151e78473b1e40ef63b313c8
cde33a8d0153f052a9c7db2ed267f1f90a8bcea7
8f900e36193a8287d3c00bfec5e8124e
921c9bb2c4dcfd37aa06725bc3339e6050307022
33c708120ca5cf695a7f5afeb636300c
ab8a14b93dfb46aa76a1f16a68a52ba86fc3684e
d0cd719d909bc70d2d5fd318c6bf625b
23821f5fcb3f398015b7b06a51736d3e57c95f27
eb7e3bac020b33145b2e5e1a9e030666
660c6a31a87a853c9af6861e15610b7b9b50de97
6103c30580897784e0c4b37c0972c7eb
f7dee398e10e5300912d16a5a0941d77d84744a4
b969a108fca616d53761455ed1fd515a
51770bb7f1aaa35a2de4390059cf1c330dc07a27
8d5269fda08a034615aea691b164580b
c67d038c4e91a405cb513e460fd661ef
4f478efe980fc6fdcde20cfbf9e41205a0bd2edf
b12507a2a8154bdbdfee2a357a60e024
c083cfc7d93486c242bc9f00ef7f6312508195c5
cb5885e91f6bb396f3bc6c078f59c288
7da05a780a7e1464e83701491ea50bbe89b46d63
968880a3b7fbaceced8fcf7aebd08637
c7f530b554183c2a43726a92637dbceca4ecaf00
338963db1b42431d8ee890aa6fe50713
6eaba720c5d0c9afa97ffbf1db005575
cb30e8a4a85f314512f42a9c7ac8b06b8b2b4a9e
a838894fcb20c0f130b768576ec16b0a
187b6e22ccf3bd6a14feaae613f03ec5cd000647
2377dfd3ab39d1d51c57a7613bd4feda
31b5d85909acda09b812b0761fb4403bef3fbf86
312919f32e02acc787be716bb5e9fee0
e1ac3a46f87f46da4bb9bb79df5fffc9a9950375
886cea20c76487d7446f3265e94ee3b6
91df7e223454d3a57255d727008c301cf19e7856
ee099825fbcbb3715fe14642e0caed4c
4f733e0abd8b7e976ee5faef7999c60d44d4d792
511b6ae9b587ee8dd907b4fd5a9ca360
dc452c85508e07d4239f6eed533b845a
e775d7f10faf673130351e2f2bef0b0650245027
19aa7b8d6310425ab6fca9083dc4b027
65b2930cd7cedcbb60e99ba15e514d719a5f17d8
839d653c6a7b4acf9e078b7ba0fd54a8
96eaae1bfebf338702c4ae427c9525c3279eac37
19d0d9195808cb773da9b167fc2d2cbb
0a3c1b96bf9d15b94ce5f33704038d44a49a3c78
405a008a8bf42385b05588323517cab7
d9b71213c0e750ff53acfced7458c1a4292a7720
b98daa7567bfe8ae021f91a1d4c2e698
c6e2297cda3e1e2b4489e9fc36207a41
d1e31f1a0c09437b8898220f654666933a3a93fa
71f44616fbfd7ca07e1e736881c76edb
643c59c8d1fd9da0b1ef337be8089e4d94bc8562
8dc9373a633b17dc04ac3bb34fbb39e6
35b389323ca75be0a8fe7cf9c86810306c819d5b
9db7c645feb5839cd8224786aa9dd766
8ac5fd6f57f1b3197c02f3a28d1c477f7655da8c
92c5f2613f3d3f75c04bd808a1663ed6
6455fced075f35fabb5fc58f26d48c869b52ee92
435621e27d742504094c7d056e967d242050a312
9fa73ca66accbf16f0fbe82b721b5b92
aa32fb3ebd74689e9dffe5f551b6018802543680
d39e602845e945867ba3136978de7a2d
c3508eab28ba3577e454f342738bd5c393db9a3c
b50b98728c40cbb4f836fc942db8760d
02596fd70944874c5fd3ef39bc37e768ef2bf393
b091f2b325f44b77560656882ef404a5bcae6600
02552e05617b36c48afce56ce3da83ac
209f2d914f6acee59fc932102513c0bb24e6e47c
a27fa0f65c0b8df2d835367bec1cfcce
b6c6a2cb55e88c54439d5cef341e4829679d89fb
aae18c2f71f881a9d6b900444ca01009
3dcc1c0cbf1411d8c16655cfafb17812929f5e7f
6b24abdd7064c50549925f4a06985cfb63b1841c
bde4872b7e07f96111eaf8d1a7464de4
a1246d8d16271801381c178fb50ccd4697e4a63c
bb1b756b03fded58023dcddc3bc928c6
7523525a05ae4d07bdddd1eea0c66a833cd9ecf5
a94d9e31413c93859eefa10b439f893b
0fdf67850de12e954a88d1e5d8bd0ef6
2fd35d110f1a15a3e0830999590e079fbacebf5e
60fa456bca8e8d9f81432eebe860839a
d6f1468f74554f2d6f125ced94134cf1f09cd938
525ee6cb588a88d69dc4e7583b07bd73
33c6d93674a486d368d6f046d3d17ba888285597
040db4f0bfebfe7cc7a437394aabb2e21360a972
6430b87f934cfc68c1d18f883e8bdb8c
7b86ccec0471818dd4fa906336982b264cd39592
78a45fb3cccb48dc0c6b14ebcfa9a332
c4eeefa45d39dd4f51ff76582415661df8278e90
6285c140131d0b3813bb1cc0ee9816ec
9ecab89a30eaa7a76ee6ee02f1f29b392d55eacc
b922824cc3bcaf3a3ec8dc54247e2a09
c57eab85c21fe2fba0456f056f79d3ab
fcbf48a884df7d2cb4eb8922f5eb46325db2f90d
06391f256a54b33c5cb0af28fa0cd6b6
2926d6c0e5463808d97ac126c59bd19e820f4e6d
974ee12230d80aaa84dc9333463b9507
1b3b71a8c10c4240e42cc215c03487f2da458cbb
ecf8f6931d4ceed2e2b52a9f7523cf29
54c243fd57cf352a42f5bf24622dd6fc76817eb0
a10122dacd68a3458f9aed7fe4f39d2a2ede70ae
a9c3048e058d5fd937c2520c7e39ba2f
5a15ce77e5f4e1b84c87cdf93966607ddf24b2cc
fc5453bda832a05cc9938bf63655a82d
1c2ea4117958f9c1cfc139319a50185923293778
4d19e8b5467ba27629ae6ce686bcde26
5e5eaa8673268f96cc989166d917fef7
1363a06a5ab8ca34337edbc428497e79fc5770ca
e81a705e224731debca177e242f98382
bf0c2a9475979ac5d4f5202463eea5d36e23a3c4
d9645fb22d72fe4a9e9934d41c4f9649
af9d861a16acc1d1eb1e99cf91bc29c47ced6792
aca3a3c4e1abb251335923007e6af346
539c85e106b36626530e1ec58825d84886fea966
5bfb155edfa3b5c3c5e31f87c7515cc7
7d45582fd947f65bf314c5a22b05ed072efa69c5
661310980ca396098b3184aa735e896109fd56ca
931faa3cdd1b82db6d6f951f82546e3c
54a549a5ef4eb7c8efb15245c25db83380ea7b37
60cdf66fdd9ab44b39fb6d075ca915632a8f0b74
8e19a9396a5d71545dbbe3326e669ac2
3a321be1052bb889ccb88233de027ef07d3cb89d
89fbd496c17180a530dd44ff43d48d1c
54cd082a243cc09281ddf783d366de64
c22220e8deccc0a25940b014f0a6c65b0ed028ec
20472247059cdc6fa598765f7d55f135
d1b1f3435212de193d34df5c98907ce0b7c7045c
23bba614857bf3bf87985004e5a8b87b
126c352d3c66fc9a8789455b59de7fd0a73bfafd
70fc9e08ca51668eea1a388e6b59491a718672cb
8bcb880197ca9e446316bcf36aa8e033
5020954510c406c72664d3d44a11087106de8be1
03f749a0dc0eb6ecfe3d917a4fdc44e9
09e83ab5a63cf4493632f3981dc43807
193d1f7b09458be5f9861eec45830cf0a929b127
b007314d8f8a405991e30e10ce341d2d
59a8595efe94945ced045284da3247ea9669367f
5ecbc00a7ae3e2bce6be65155c162b46
884cf18e0bc8291df886e150e70d2608174a41e3
c5f419a0c750d6893a5440f866592c03
857d44954d4d06b331326761b3d23c720b647925
e11c65c4f0db09932a5204cb3416cb1d
a8203cf02c9061268119db6eb95416de12bcf7f6
570a43e017dd016428dd36104ae18d65377c8bbd
3417a82bc7a9b96ae654fc2ac02bfbe1
4131f215e1dde207c29b6f8a4fc5a2a2a3c86748
c1c77faa5b2dabdd7e261fb7b99f1c96
e665aca6ca9f6a51fd70806b1b18d6f10b351f04
d891bdebf48bb05345490a287a53174c
818e8aee7c90fbdd1c1e0d3d61cf6ed986579e8d
7cd93ded30bd5808f21a39ff53a57282
b2379252009fe155c0c9749168ba147ac3b609d9
e61c97b7e8d1c6986fd47498735606bd
64ab85a34b6bc70b6dd3596864e4e55e60c3370d
fdcf6b2689bfe9ff3faa2b7407d76e44
e8c8ea56ac313f4d30c39b69249953f9ecd65eec
48e9eaaa667d88fd0c4968999ff88c88
9ff05a7ba106fb155fcd990a461d80b8
57af3b554d3f069cc25a40cfc82d28ba
6ca5d6e50a174ae338a65b6ed796acb81e555dd5
b354924deda9b883f10f64c587674dd9
9cc8176ca64a8c7188f0be62f7c51e645ec4cd57
8f3afd781367fd86653556dded0919a1
e54bd8f2a9e552b4c3a967782c21aaf1
a23421fc793662f567e0545a6a66f59a
a0afb87d2a80814a2c6d2733e4132cf22aebd805
2a82aecdb89bfd05a7e5dfa0ceeffca1
1b12ab5ecbf00df5fdf2dcab2e5310cbb5c6bc2f
ac3403326d028ddffb97787a877370f3
393a6edf40a51a918eabb80996ca957364daaf8c
e05533d2673a64f8e091e16b63631a91
ca5a3cb606615a7490675c21fe2c77e30f82d75a
db685832a5910d4ca5093495da2d3e166d4a4d84
bcb404d63743b8396d97c4533f5b9cd9
1dc7454e8bc6a95c7b4acfc0d3d8dc00e9c51860
3e56998b91320dcbb1d913a90662822a
bb85a6226ab5a1be1a44fa8a607b6f471a9d57de
eeaef506cc96671f9ee7edcacc51f1a2
c8283640a2036715a1e985a8b7c02944
07eaa986030fc4a79c0e31776b043c95c284950c
16492bf57cdc90885496fd104b5f239a
43c267bb309a090ce06a2958614255a844f556dd
5bd139aaa55ff4ec55cea143a5342f8d
b320494a8e4817cc80a6c12bf1f113043e0efe2f
a5213dbee75136fdf072598f51891895b00c342f
f8cb11df66001d61149145655de1dd404f8c297f
39884ef69a3842ae07342802a45cda8e
4287942af05356637d150f27d27853c8ee1188cd
c9cfe8c73014dc11fab730e93d6eaac5
cb286472faf6fe696f1ac9b003148318
15c2c1b3f40bef7ba6d18890f13a18328c9abfc9
f8e19e5d919aa454350297de681fa303
b32ebd1c33332700eebf693cda69e2d999fc73ed
b68e8fb91353f4d503ca5484eb31d3e3
5fbf02b4016b193972bd7b1c713d8871a4f6b478
14c83ee48fe155a89b8495f016d03ac7
dddb75ec495bd82cedda2239afd1849afe802d0d
37fe4c8623b42eb4c55a0e7b5cb021fa67b232aa
30e59003bc794e3b75272474f8b96149
a687f07d1be3e1cf5c88e9b5f87ae727f5a3471e
fe95962036e36ca3ec3b1210cb030d9d
c21bf1f8a5d1dff58f4015b5ae8b3fb6111737c6
7fc453436668a3fb89a39a60b2918454
598aa49318bca4c47b8749b8170ef9eecbac0cb4
64a606023b9092205f75e18613d64894
790ca0fcff334a71ff12f1d60f593fcecd145bfe
bd58cbbc8c75a0a246b3841864c526eebfe03d55
877d6de64afaec67b05b67ff50224aa8
edd49f58df1e463a4d7a6074eb90fd4e6b0d84ca
0604d800a769c99943d63b23bd5ef03b
c656aba966b711689bba7e7464f946bf
5564aa086effcb1c593f089bb7d8d940dabc7097
205903a80c920d935b261e2b162194ec
697a5818f5fa5b6377221f1a3051d629477d9f48
32efcca7eeb8f687be1e14562e43b220
85a3409674effac75b8286ccd5a0aa79245f3b27
dedae96de01dc17e30d74629708d2659
9a247afa8bd253b7a9d13d1e8e43531dfb39b7ab
4c857f4bb889ec5c67b40205a559a143
e672f9817a06955612db3f0ec952dd4fd4d1b8e0
49e20e63c2a000c4978489189d4830b8
6116f222b1d493954e04d12ae9dfdd92a5a66280
f83f0077929c0cf4e6b7967b4f3f87f5
c3724f2283f33ef696d27e9dafa1ce06b7b7936a
8729ceff5aa3794ff81baa057a78ef9c
e286628e27fbb3b371ec55cffe600d37b4e04d21
7d8688bbf8d4505ca6ee309aaf6d54d4
886fce0f7d1f865917a2d30fe62fc7d401797c9f
b56f776a3d8e1a557fde9d57cbbae256
c17d40eb142f25ba46def399d10f211b82b1fbb9
43305b2de185e559a6b5e51f252f1582
5669c3ea2bdc366a8b605a2c2b858440576ec4c6
d25026694ee7dbcbcb5fcda2b303acd0
6c9abcc26507d4fde61a318ab3ae0d1097a2399c
25da557d7347bd6d1b2c9bc7208dbf1ee7fa7447
926172ccc2ecb74b90dd989d87f3e821
ed3fe736828a8881bf2323a7f5d48d3c9ebab3b8
14d62af1c2a130517ef3a1e090ea4ffe
22377538874f5b7deff94e446034f035
21fbe86963d5eeb4d3d556c83b063ef0600fd38c
fd0b506eb48e0a9fd63d78e93ea70ee9
fed8c25bfb884ca1c565c2daf77c8de2968d36bc
ccbf9b723885b3769376bee31d29e462
0e1e422ddbf017c19444a9ec368ac5298a6500b4
212044a1c7fd9ff896a0857eb29f959f
46bf1b3cb63479f8b013a063f13f2bf7d94c697b
1c1fa35883e394a8ba32558b556f4f45
aea0ca8a0e1c2469976dd6158c93d25f8d94ea5b
e15750b6737327c7b472f147c6a44a1a
e053e397837c453e0b1510da1a876a706a41157e
e2011fb8b4b5347a559dc0c3996feedf
4c3a335e72e60e2dd39cb4c85773340590b3af3a
387224332a810eecdb4dbc84d5c7a779
4d0bef0c68f58067efa0537f3e7d485f29f6b72a
dbd85098d4360d2818e475c7c3136c80
8663384315ee6bf5bdd495d6325857c82f73542b
43f47430cb5e95ea85c312b9c6268b6a
006d0b81a944ab28c9ca09ef3a8852ae2762e9e6
a29b41e6512a8342bfe4f83dd84f0373
22682915df2ef22ae970eee243ba5b613ecef93a
5e646627947cd7c2c1bb73f66ed567b5c38a0779
32976f9080d795c2bdf19014ca709f7a
9872fa1855844dbf71029ebbc3476b6b
1693b3570a5638cc835b1701f829312793e7d65f
813b3e420281b150a74da2a8da58863e
a80356afd2353fa84a627fc6ab1f80f178fb4d82
db12b76dd005e190d6eea6689b6cd6cb
6f8342b5f98b2ce4612284e68ba472e688e19282
2b3924af79ac0dca41f2b8fe14b78e77
4367e2b929553a2eb1cc2e90047904d9423403e4
fedcded988eb25df8c3ba67c3ea12921
c084f7a8cbe924ee9a5b4c108a8c3236d7f54d63
c465bf98269fea442507eb7ffcdc8fdc
fea7ab7263c5e18d0c9b0e9fdcf3cc5cfecb81a6
ecc6be6ca1063bedb7aff9e6623f0209
9d4567f763c5b1942de2e643bf9ddbf0d8bb3ecd
84d5ed90aa0b793afa13dcf4e141a74eb6286931
e862ccd2fc908c46ac9c159abd68a2ee
17230252037684cec7bc1b6ca23cd5e9b868e7af
fd0e6997805dc2c6cb721d6db522a360
d1c7538681645d9cbac84e1b87580c46f4bd05d5
43b95e6fd054e63b7a4979bfa0a4dc39
69f63bb85dc64368914951238fe7756ca93584f4
89e19fc6a15dd7320bd93c9bc7ef46ec
f1fb22b9290a407010e4a4d4cd2bb6fd9059ef12
6d3b848b5bc4ca83639f4f7ed3e4eb75efb086f2
14554f6119ebd375feb4ef0089a670a0
ceff4a47b314db46426e71bc4f6384df
3aa6f9e8c2baeb9e8c9ac423fc43dcd668d53cb9
009e937813e1f28c009fe946fa159f7f
37bccf48fec392ca17aee29c2a4ec7ebed01c9d3
01208ea766986f64ffbea75ae17de1d1
1cfa00625aa1a84363c295f10733f8c169312322
43cac5a6bc531be9826ecdacf044b765
60d7bc08efcdd1c015efaf95b02f5c230ba7c994
22c9de590562f6bdbfa248dfbba20f7c
ee04a0a19bcdf50373b29a30fda09c2633504ef3
8ecc484257ae8468af43388093dfa175
081f49863a824d991140d3da72390c2014a8923f
b90503551fda8f1ab95a76e7f8149605
a67fea524f10e006df58cfa208fb52b65d113717
18fd3ec88030c5ffb37aeb54f1bd1322
1a7dcaea093f803bd55575cc997e1b5a36ab873d
4c6b2bfeb10e1a26b94187e14d78fe3c
80ba5462d48826b1bdc025f068369bde
3557e81ab77e6dab38a5e6b61efcf82d74aa0212
7dd5dec0d4c609cd6f2f31415bdab5b0
161c64965d28069b0c13bf0cc4c5fa86
aef263b81911242398880a753db31d5b
2b87f151d81277bdd340e675662f3c0c3d4c67e2
2c7cdfda9df2afaa354bfbc6df896c56
c16f6aa656a9c446c24bbd92bda5f4631e961688
3ea11816df22748f89fb06651ed22d86
764de6ce910cf81789be49bf987a165987f4df11
09796083d7e09c46854aa898897ac9ff
d1a891ec68ae09c5eebd0cfa0aaee8fd4917972b
0c0b3c171453307cf971c44b1d344a21
370f934b12fb7eb46d9446190dae8fb3ffefa812
afc81fa5faabfc3e1af52270179e87d7
1447127a94c08fa49a8645b40a594494
bf239d408a02e25c47e460820fa45c5320d9dea4
eb64dd308fb6f52256cfcb7782bbe288
16ef91a12eb92a08b607fe2c0db6250bfed7d636
9b1f15a8d283ecea3c8c05cacd7504a5
053ce922b80eb656bdc268745f93c614c798249d
ff9755ff6f2ff1d937d3008cd0e550aa
c3821758f05d34e902d1f273b0ebd1196928aa10
34fd4426009af0ec2cdc449898b619ea
8cb1371b941da0380ae20b72d0b3ed8af44645ab
72153d6d33f2e56ca06e9d7d5aab9e32
84cabf3807d56e395dbcd0705a1a087158bfe13e
16c3236702a0d754bd455a1f13ee3b7a
7170fe14ab5242b23704256e7fe081ba3b339f48
1d348b23ec4e3356fc98cae5ce2932fc
bcbdab8e2ec0408e7e8aca7b8f518270
b02036304c3b92c492ea30d3c5565eda6d3a6333
a2eba293d9594e6e8bdd0dc3f514594f
e830f46662c88392f9854441b1bec23ab2e89d11
cc0dc17a6b2af49b02304715cefff6eb
7651 F20101221_AAAEAA roy_v_Page_14thm.jpg
9761d2e535f47c19144038213472ffba
8eace225d306c700811dce13eea4aefda1f27d4e
300ddb7a2dc807a3ae96b97dc2f34911
5b81fd07d77ed1c76f1080e9030f32e5743cb196
835187ab785f43c463b3d70ee3dc7943
4d13d6402da1cf9932bfd061f5588bddcfc4514a
35463 F20101221_AAAEAB roy_v_Page_15.QC.jpg
f78967247bd64b02bcc059421e1593fa
0e4f57c1b4c3d90e68452d005c5c78c7382b5fc9
9be9a4f94dbf2e868083961a5320f914
6f39804986d7b2dc0d3031ff5786d120aa92f8ea
34354 F20101221_AAAEAC roy_v_Page_16.QC.jpg
6543cda96d7ce9964d842211e9828d09
44a5a50da123d3fb4dd2a91ffcc8aa9e
5133184fffb1c41c86136e257945b5008e22b80a
a1bbdb39e37d4b5e3c3035e25aa281b5
6c0433bf6d4a977c727a5a44a75c6fe518d23e7f
85b09989fedebdc58caba6b6258bd08e
722811c3b1a0059d98a1df807e0e979e
fa0704a09c47358f9c4fe7e19a916c4ffc4bb29b
51c42117c5de714be3eb43c48b78940f
919 F20101221_AAAEAE roy_v_Page_17thm.jpg
39b57d934db3ab0712bbfe2fc032520b8bcde4b9
dec53ec4d905a33c9d53d01ed9948696
dcfafbd232665344f744f98ddaa7fb340655947f
0027b938008babbd49e5386df56b1abee2eb7743
21233 F20101221_AAAEAF roy_v_Page_18.QC.jpg
77da661ba46e2b0eea2f735bfb7dd9140a44cf43
e1ed5a905c13f7333c01db7b850a6dc7
df21bcf6bdc5301894becf99ea9101e940f529c6
6d0bf4e691456051a7821df138202eed
8e91d97406e281d077c5bf2071944c769a597f52
5493 F20101221_AAAEAG roy_v_Page_18thm.jpg
11dd85506c648935a565d4daa27ba250
313a4a0943aaf9c542273ff8c7b47584e6362160
80dbf3ba42192f07bdb9c42aaed8c862
af8171b1b5c16dbbed24526bed10ebf12b333eb7
6160d7b65ec2c0aa0bf96250725ff916
e719066901a5453425d5ed9da76d3095eec1935d
09fbc7a9445faf5e121a959735071d35
073e7061f77780dc7a4ce92b8a14b9c88e867e47
28103 F20101221_AAAEAH roy_v_Page_19.QC.jpg
d0d601e27c1d128f9ba94ed78f95b5c2
96babb032af9339a5544568d96d27d24a9d5c822
23cf69989ed230fcc8b8fe46b3d61b74
735b28491337802008dd9c4cc1168b9d7f047407
9eb549577f50a85ea0fe9e27e061dd4a
cf10fa68c42e1417580a3d21407b62ae27225f3f
799ff6d2f33d018be64de2ff8c4028c4
252c2d705c6e3d7c610a1298f681618a1c3682fc
27427 F20101221_AAAEAI roy_v_Page_20.QC.jpg
224e545a24d9e81a3677a16e6b76cfe6
949e6dba906e387bdd85f2be7004475dcc5a5e2e
b67fc7da170df44fb3d13cda87129229
abb18a8b9fb5f553eedf927948e0e1193f4c1e12
31ca251c695d8a293cf27724c3d49560
9e65f657dfb016dba3ba16b04a24deb6dd0a39ba
7070 F20101221_AAAEAJ roy_v_Page_20thm.jpg
6cbe319f5028fa7d78a587ff8d3ce05f
62cd8f2685943b2484fe4274252779a8cb6ce48e
118f25a924d836769da82e221b054c4d
200ca94e796a7f0e9d7aab47fa5978502fbd04ec
2b860cc06bb861956101f9a0bda30862
2ba462e706bfa08333b0870231612dbfe2940b5c
32049 F20101221_AAAEAK roy_v_Page_21.QC.jpg
81721662ee5621aa0df12c628144d4ca
b8e3db55ccf849dd940df059b7c080602f5ce278
3dd877432a7997e97869eb215f68536041f41ffe
5ee5ec117475843cd51d44f69e3690e9
3b9974854539dd0a13a86a4f6f8aab7a9a379334
8495 F20101221_AAAEAL roy_v_Page_22thm.jpg
79f5ef3070ea78df4b2e601165dd5c4f
885868c54f25fac73034a8c893a102d3284a1092
fcefd5096da0d0bae7ea4c6aa3511823
4f601d75323d5c8f9f007b32480ab1107f48b62e
4bb32cf1f35576b8d93ec402d11c9d45
23186 F20101221_AAAEBA roy_v_Page_31.QC.jpg
fb1fe74e50e25178e8094d076d06e183
4aac7b86945ff1c5c79ffd692cd97b7490edda6c
19620 F20101221_AAAEAM roy_v_Page_23.QC.jpg
687296b32bdff944e27fce6810657435
d2165e0d242b4f1434075a2bcf1d1182e3cc0410
df27905884d53e2595e28ce96082a154
86cc7ca30f86d33994d69dd831fc0df5e6986efc
82c990f863fc9924ef43185cd070209c
15f851813822e363437e19e2588756883736ae5a
4678 F20101221_AAAEAN roy_v_Page_23thm.jpg
aaa7718679a2d451d22b8580d6e645af
9f68186c22c0e35ccd9f6533d0020f8863461c10
bb3bb67713263481fb82d9592004d56e
6d47a85d4a4b07154715638bbd74485610289096
fb6e1c07dc0780761dfa06038dcdbbde
6426 F20101221_AAAEBB roy_v_Page_31thm.jpg
7c74aa06eabb9d5ebe3524c36bf1c9ba94374f15
29963 F20101221_AAAEAO roy_v_Page_24.QC.jpg
964356d84349422f04bb4ce24e0d3fe1
ce092c8a7be2521f924c03cfeca65a50484b3a48
0b92f2d17e0cd6c483b0df1c4d1aa4af
f8a939b0960b0b9341b3e7846ba64cc70d43a500
2044cd67c0123406e2b747444bb6cdfa
a5dc585615f3271f5ce3a46fc7dab38708556b12
6270 F20101221_AAAEBC roy_v_Page_32thm.jpg
2d7f5c44d161d342a087fd77457e91ce
31343 F20101221_AAAEAP roy_v_Page_25.QC.jpg
2a8a7f37a108b553157c53b2a290d84d
e4d5a391d48f62bc8161fffa930d6a3b67450b31
6847ea326e1903e997c16b241f017c2a
e335ef593ce0c76febc83a07aa7e0c05ddb15315
17378 F20101221_AAAEBD roy_v_Page_33.QC.jpg
873d7599a5173a8d7ca6c420bfeb6cab
27e6868b42c01401a87a542d2bd7210fe97d7da8
7558 F20101221_AAAEAQ roy_v_Page_25thm.jpg
eec1eed95389f6c3a119c55bcc9c9a4a
5d36b885017b555c362e51c552953231df6f1d6e
d6b09a727abdf9f0b4fc0c8c613b0653
ce86bd38c060e9bdbb70e6ec9eaeffddc266467f
9bfe23fb978de7d3e720b265ba0a3271
3159b7a15ec7360cc3266d24b66a2a6b95e994d0
4897 F20101221_AAAEBE roy_v_Page_33thm.jpg
41286175ef00c5b15c8df17158ec1c27
728f4b638851a712b2f6a8de854e465f1194b105
32032 F20101221_AAAEAR roy_v_Page_26.QC.jpg
b7e3e5f5477b0dc21da26bc8bc472a4c
ae143097cb608b00a20765ca396350059942ae62
d077f82ecbf6a5b93ce9b46832355e48
5aaa03c46ab800792d639e7dd0154acba73e5ee0
5dac65b67a98e58fa9d1a7304d72403f
011f6296e351385eca376d0259563bb5e82f55b7
5168 F20101221_AAAEBF roy_v_Page_34thm.jpg
3a6724cb8f616d149234f2ac2070a685
46487a032940a3a1118a58e2695bb397ed9351bc
7861 F20101221_AAAEAS roy_v_Page_26thm.jpg
9cda4e691f20153b39e78bcb7f74411a
23606c41eac30ac64e902ec678f17abb
9ca7bc7d33cf4c1df59661558cf1c80bffbd133d
862904fb83bc845b09bd0059a464f1ec
471be599a127895412caea41ccf19a968f1a4d98
32507 F20101221_AAAEBG roy_v_Page_35.QC.jpg
1681fc055ef5bf0a3b0092d06153f32a
46909894fe92c5b4870f62d6b46948ee97e02751
e339ae8ea363a8700207e728d1aeb17e
28154 F20101221_AAAEAT roy_v_Page_27.QC.jpg
d25826f3039c01b0487dc7c96f36db45
6f1d3dac39c8bca4727de2ec037030fcbb519873
30a1776d6bfdd165409e96d919ce66e8
2b4e2e304444ea0161df0082041397b146da63eb
a116b5d6ef022b3cb9f6013ba356db66
db222d3b3ef90872011a6166cab385792ce02d34
8011 F20101221_AAAEBH roy_v_Page_35thm.jpg
c349ec940b56ff437c230f6bddf7da6a
3d3a49e85755a743f299edf1973d32b17767ccfd
fba33857be850e80f7aebc8b8e5575f6
5abd402609a604f60203ba06da914b07b33cfbe7
7002 F20101221_AAAEAU roy_v_Page_27thm.jpg
bc5a7f37bf003e7103343cd36798d5b1
215d0cd11b13344ab291befeb2e4b154a68e4400
e24fc973093993e43ea3d684c3ef1f05
4f963312fdd1c4975f64eacff5869ba2e11d8875
90ff8e3dc8838613e126cc8fa8efea9a
035ebebdc2f87ef1e50bebb0226026c6a79b3730
31050 F20101221_AAAEBI roy_v_Page_36.QC.jpg
f82fd8e7cbb7e593c2d9e18cbb017d1a
0a7069e9bd76e35e4cafda39272e920c667140a9
6a15aaa7a927fa2f6887735f67897267
8664db4237b7fbcb824e47eee1e9aaed4415b8b6
28130 F20101221_AAAEAV roy_v_Page_28.QC.jpg
ba9fa000051f092fc019f72602e389b9
f7557f36d78278b92a808f86e403c4af76290abb
a37d47207f3830c9d5d926eb609f08ff
626d2201380e0dd353893c966273a164a531f48f
7477 F20101221_AAAEBJ roy_v_Page_36thm.jpg
6667bd220dfc8a5e6fcc3c63483bdcf3
295b2080e23d119caac7cc5288263e13c84e7370
9a511c0633147ca34abc882ffe7c2423
8d7816e5b2d87e75eaeaf6c60eed275053d54b69
7144 F20101221_AAAEAW roy_v_Page_28thm.jpg
68d66b3b31861c0c518c9a5f057f1b0f
9d6067707c0bd84e21d421626412a6f9
7db7ce39bec2b1e9a71d29bb89fc6a779f778646
F20101221_AAAEBK roy_v_Page_37thm.jpg
fec54fe32d4b5ea5897db304fe44636c
1848b76c284a202ba0397e35a7a2cd54c0882c2c
5bc439bd130a47aa1df0e2b81811a79d
d13ce47f6a97616a373c86fed5abdba798aacb36
31669 F20101221_AAAEAX roy_v_Page_29.QC.jpg
fe98a80402ab2c286ea679a37315e1f3
6542f820be473dfc043479c36a96809f4cbc8537
ecd550b1c9059691f55bb13e9424b738
05a1f626179dd40e210ae1ef4064a7128ea874d4
6412 F20101221_AAAECA roy_v_Page_45thm.jpg
2b6675c208d71fed6086727078c17bd6
9c5d45f1c07c4a296f8f6ce87cc16e0b9cfb6411
24117 F20101221_AAAEBL roy_v_Page_38.QC.jpg
a6cd37c6f6d0da677490f70df3b08f55
4d3163313299176a67aa8a416913a575f63a7d78
7937 F20101221_AAAEAY roy_v_Page_29thm.jpg
89a6f58ea6aa21102de31166fa9846e0
e6831bc63b8d33203f9371cbcd894eba864c7185
25862 F20101221_AAAECB roy_v_Page_46.QC.jpg
bf4ddd5e38826e6d8188e8ea7e9cc493
0ebc78504d8cf9f7ceb21857d3efbaa0de07e2f3
6499 F20101221_AAAEBM roy_v_Page_38thm.jpg
a3203cf08fd8032285a7742e2f315887
b6389b264f0f9e0cae7db7e9327719ed7713c0f5
6ff513c50c567e7b02dda54e2931eb74
d1435ac22dba56326d47107e3a2f94a464016257
2e4a5cbb0606feb6feb32376a4683420
5df209b74cc49ee7676115d1f980734e619c3a81
29918 F20101221_AAAEBN roy_v_Page_39.QC.jpg
53d30beb5b2f41ab71c0c1aa9b6bca21
30c68a92e97aabcc6afae75df938c7f97148e766
20098 F20101221_AAAEAZ roy_v_Page_30.QC.jpg
d8b47208c5a66ce85663ed1649839025
8cc7c2b7cf5ebe4dbb9aa23fa657acc62e5b4db9
65b2fd053da1041197ca013be1e0ed76
01ac0dfab1882b6b7223507f2c74029a54981efb
6474 F20101221_AAAECC roy_v_Page_46thm.jpg
bc07c127e0ab4881a0d46b8c93051bef
f35d6c66af74e4a99f446b179f98b94fd70d3cc4
7120 F20101221_AAAEBO roy_v_Page_39thm.jpg
a88225d56da84eb2b78cd6a634a90a4e
911e20a3d78ec6566230d3de91662a05302aafba
8908085ed35739d59d10eb920a119bab7d34163f
202cc8c9fcef6a23d3c9579e541b94ed276b46fd
33275 F20101221_AAAECD roy_v_Page_47.QC.jpg
9f6a3f2d3547bf6a9d79859eab9c033a
7f6480c521e0c4ce1fbf7f8a5dfe7a79505e3c51
26561 F20101221_AAAEBP roy_v_Page_40.QC.jpg
98f354f504897e234e356cb3664214d8
b127e2d1d551b20de26ab5fc49636a43027f0a9d
9f5ffe906ae3ce477192472c2f345002
4873ebc98ba70fa6903f107b4a90549af7cd0477
7847 F20101221_AAAECE roy_v_Page_47thm.jpg
f573447850962d8df9cb02bc99899bb2
322f17e763007391bcedbb115815d1004aff9142
6813 F20101221_AAAEBQ roy_v_Page_40thm.jpg
8358136961de3e2652e585324613a986
931e7d9c37c6aa1f4ee255c8738b3bb85c939b25
65368a8ea645acb23851c74ef88869dc
6fcc0d7b38c570f00a42235e8a8a11c3d7c88b37
cace2b6eb3a6a150b9a315e2f309f574
e029aa0be20de0ed49f5782fe28abb96fe0b8813
27087 F20101221_AAAECF roy_v_Page_48.QC.jpg
a42b58a3136dd0f71b0e3085792ec34134609e61
24979 F20101221_AAAEBR roy_v_Page_41.QC.jpg
e33dbf4944eb0300f899f6bfeec0aa27
05ffae42e2f248d83b1b179e141516aa28528af0
807d3da82872a78f48034715b3cc3347a8dd74d5
4c912441c670b56ba1a14332206e36966087fb3c
6915 F20101221_AAAECG roy_v_Page_48thm.jpg
506d1df45e2d551dcb49af2c8d84444a
72ef490b12128ca06fd29b3f6898d737ceb0eee6
22b9e21c8bc23c9f45293956b094a484
502ce46b8de7aa4244ebd083757cde9ec0a7a7c1
6199 F20101221_AAAEBS roy_v_Page_41thm.jpg
f7d34d8602dac134e8a9cb499fbc39f4
eb132bf915c787f235c4cf4aa47c3e59e197f0d5
4871c07bcb9584cfaf474a848baa1c65
d19d7fafb3c72bd2b4e84a3ffff443a658e98b08
00786d4a1ae0a9e899abdc39b5971d37
51587fc37da13eb1a5588b0633dc7f2f6bb5cc3d
33205 F20101221_AAAECH roy_v_Page_49.QC.jpg
d46018a0362a6b104d14dfac87466f7c
f3a92dc57329ae696cf0a54df9670e30d236cf46
9753441aef6634d784f717dfb4c970a6
362b495979c3ab258ba0e7255e78989410ff094c
12703 F20101221_AAAEBT roy_v_Page_42.QC.jpg
74d7147713c0fb530e2607222cdb12a5
7e35847a7d122beb1c5c926b9cef5f4ca3cae9c8
a12aefce94eb80deaf71508ac6211566
153f3db0920c03282decdf54806f07b5fa89fb23
1a901fc7f25cc9d83b0c8edf249117c81e410919
7799 F20101221_AAAECI roy_v_Page_49thm.jpg
1f0f83494bdffbf1a5b2ee404cc667d4
28f3c2a8160d5156371dc878f7f01255
e22432912d07464482a1bd3cd664249702790bff
3335 F20101221_AAAEBU roy_v_Page_42thm.jpg
c0ff1d5e3ebdbea90e2335cc557f6800
28b0bfac1db6e5fac0d3da6d2d363ce504c2d8a8
b4d5a445669b02ed1e0f530689aa10da
e3038e77361223b40201b1579408bce47e960fd3
32547 F20101221_AAAECJ roy_v_Page_50.QC.jpg
e1fac6106bd9884800e01b1a748ee594
2b69f01da36f6bf00af294d3200590767580ebc8
91af8059d928f9a0f64d3a0e73bf403a
33c3d74eeab6df9c90544dcf08a54cd3d81e28d9
28980 F20101221_AAAEBV roy_v_Page_43.QC.jpg
1a137e78af6b8beaf523b78911d76ac1
69758622e73597ccb09e0d00aa916455529938e0
91f2eaab0777cbecc23ec94d31acc0ff
bbee0f3170b5f12ca868d505f68561feee5cb421
14319 F20101221_AAAECK roy_v_Page_51.QC.jpg
f4a012ee8f7e99e3d00165c76b6e1e82
5a41b9ea4e7c6f4eff1dab03d6f38e41c44a2d62
7010 F20101221_AAAEBW roy_v_Page_43thm.jpg
8fd8784354c33235a5d288262c24ab26
11165c6bf721be4f195517de659a7e0b
7a2b9b2cce31d72791146e1b728e7ab5c22b674f
4061 F20101221_AAAECL roy_v_Page_51thm.jpg
16dbcc9d0950795e26b4603054ace05b
fcb65d323d147015fd0d6d07a1521d33c3322d89
8794bdc54ee069064bba17f399a8fecb46d7d2bc
29483 F20101221_AAAEBX roy_v_Page_44.QC.jpg
dcb135e286cc3109e16d264df44b7fc8
9e58b16340a472b1f9128b19363c25abf723591d
129d8d1790a62685f9fc1b62d76480cf
fb9862e978e0cbc9af9bb852106bb13e28f92fee
28647 F20101221_AAAEDA roy_v_Page_59.QC.jpg
a06ac9143c22fdc1e2703ee72376395f
259f24ab74251351b2f858899692fa2703e882ef
24989 F20101221_AAAECM roy_v_Page_52.QC.jpg
9e9268a51483552848818aa948f3286a
cef66ca3b7435c714e19bb67e2291723e687f733
8351eb33da5657ac53b493129973faaf
518f3342481ce6c4224268a5e1ed9e4e26f5321f
7394 F20101221_AAAEBY roy_v_Page_44thm.jpg
c1ae9b410f7d4f8939f336cc828783c1
2c64c41f00465b3d68784f0f0df6d2344f508d93
fca73e0e65d356e5074faf8640b8ee7b
e5e3ee479817733959f2b8507aec32ff6dc2ab66
7074 F20101221_AAAEDB roy_v_Page_59thm.jpg
da80285d612e55736a3e4f28c23ac97d
fbc95976621fb831aa0d908ecd7116b1056b9446
6595 F20101221_AAAECN roy_v_Page_52thm.jpg
2e2d85102795ebb90b71d37d89e32864
3dcbb4180ab7267529135ebe323ec010
4024032084b80c4a7e8e2db5842385d3c82cf09a
24589 F20101221_AAAEBZ roy_v_Page_45.QC.jpg
f06e2dd5d621251684b643e411961595
b7b9508430852262215b020ddf082bd3a67ef509
af29287fc736fee19c8ce05d835125b1
ba7807461f34a981d5c0ffe74b7212680eb7d1a3
6391 F20101221_AAAEDC roy_v_Page_60thm.jpg
3df352611f05bc42c3c2263416b93f70
0715ec2e79d565a1cb5d3c2b948a92dcc1e6a1ed
24956 F20101221_AAAECO roy_v_Page_53.QC.jpg
0c7dceef40bce324e86d4589734e2d11
a80466f433eb4e24df9c0d4e5f707b61f2e97e82
667f5b48ab04bb47885a5c4e918c88b9190be964
f432802c773217123c25c0a265322d55
1d064d7aec667881476b664e31d588e591ff93a2
6067 F20101221_AAAECP roy_v_Page_53thm.jpg
458736897b2e4074b4950382ae53b794
ae206cd1c1f1d6f5ef2493614701420f2f7dafdd
87245dc8ea97d1e7981a2f8b2f6a5eb1
579022bde1b122fce6ff7c99cc7c361137543326
8ed97f971a3174b34df35760225286a7
442419864263ca3c9a24555fdc35f4bc9a0caafc
28449 F20101221_AAAEDD roy_v_Page_61.QC.jpg
72a784d7ff174312ae23812f09f8fc641e9f0573
24608 F20101221_AAAECQ roy_v_Page_54.QC.jpg
af4fda1bc821aea563c25664c39f51ff
19631b02c2ccf60a3160cfa49285bd9881226b48
6eccffbac5709f0ce03a1b5ff52e9544
c3823bb83790d076aa5a6f2ff1001ee7265316cf
fec04435eaa1c2fcb8a5e8acb6c15d57
574731715c5be44de1835a82a31d26aac6e1a7d0
6963 F20101221_AAAEDE roy_v_Page_61thm.jpg
c8225ecd1acd8669fccca27f2ccceb2b
ef5af0ce83fea977814a2b55c9f331b4248225f7
6933 F20101221_AAAECR roy_v_Page_54thm.jpg
0b245f8b2ac7e8db31e2598dc5dc0d17
30acf315f13d07343606ab24ee2b6f384c742c35
90cda251ed0530a7619fb060d2b07d91
18973 F20101221_AAAEDF roy_v_Page_62.QC.jpg
c6476fbbab271fee401945232e1ebf67
b2d5088a4d55b149ff2e1d1a7edf8d56aa76261b
fb0cdc0a707d4dfda8f89c54ce0c8ec9
9445fc63b1bc85f096e30fe7bb66f24d2f3ce669
29650 F20101221_AAAECS roy_v_Page_55.QC.jpg
9517771642227d963205f33d02112a4c
e7d0173eba8246aeec66efc4d36397874eb11c0b
4868f4494cabeeb339f03aaeae3c8a8d
2e3581789cc940fa620bd5f20bda6438d120f38a
07bff11aec4877f2b8ca2aff00667e1bb5ccb700
5793 F20101221_AAAEDG roy_v_Page_62thm.jpg
ec53b4d54147baac989e86b0c00ab559
2036a473f176902fc9027d63468a73a128fc66c4
4c4706b04297e0c63714986fefe8941631077968
7605 F20101221_AAAECT roy_v_Page_55thm.jpg
7f88cb0a97a799df94d63f2eef6c26ab
64fecd68c42f4aa43803b4ec4886ca236a298510
24cc916edae8b7bc36bc4a461fe331d2
e47f5c8ea221d1d34dc48628908f6fe73990d36b
bea7112070a4b937e858b1d5661427f2
307cfa9d3e2e43b88d45e5e7b2a4d25bc600c514
23660 F20101221_AAAEDH roy_v_Page_63.QC.jpg
63097c06bd1638a05a88bebdc30e58845e797c57
550ac11091a036bba687372b10dae3c5
79dbe0a464bca2b9bc52cfdb9a6696c92501157a
21439 F20101221_AAAECU roy_v_Page_56.QC.jpg
7bc75f88b69f3a7d772ee43412d6e2bbea8dc6aa
cb6885ba8bfeb7ddb1c7ac5643a85b26
146b4bd3fe83c5f5c7df74053012c27bc5001573
9e6aa38abeaf6ba3a3a9c5ccc47a1d65
45cb725046505d9bcaa777364dce347aa1da69ce
6459 F20101221_AAAEDI roy_v_Page_63thm.jpg
e3e6debdb50a2be10ffcf2ba9c0bd7d6
29ca382b2758a9be9741781f4094c498130738e8
874a2e4fbfeeb812b6ecfa37d0fddd4d
ee43103db1c2ec763778e6d642ce5950e833218a
5932 F20101221_AAAECV roy_v_Page_56thm.jpg
e9f5572122321eacdfd3156b66133355
fcaf78a9d3f5c83fea6fca44ee599804c2a3a24c
a0e97c41b78575ed598c0b3e274ab1c9
ff98b3b3ec34bfac28c38293eaeea5f7a1393e27
18382 F20101221_AAAEDJ roy_v_Page_64.QC.jpg
3711b16e5006c7f8324060431ce788e4
6e8fd3bc0faa348393a89f082deef8fda0e9418c
7357f1072e00d7d9907378499f6e74fa
66d0f237c52bc339d682a4a530ca06f1ae1d2643
24027 F20101221_AAAECW roy_v_Page_57.QC.jpg
131d7fa0a610fd1c484710d5ac1ea81b
5cc069a5a847cf2eb3e7364e33b73fa2b35f27c9
fdd209a0488119a46235a8fd56056bfccfe23d6d
5295 F20101221_AAAEDK roy_v_Page_64thm.jpg
61b7a5c5b8574c7ee491e07723805451
7a2aaceafe357daececbb64eff2269b96e10845a
5cbaf1a5c102112ec143f5c802a289c0
a3f24c4e51655dc0f490f2e54b52d0fd07a4a0ea
6235 F20101221_AAAECX roy_v_Page_57thm.jpg
9b948753867acce257828940f5e151a5
033fa4a4cc9cb05c7bece79ff71610957c34dc8f
9a57d542e9786c7ec139bedcb77b934e
dcb8edeb86bc3ab5e9003c87d40725b274f51d97
7194 F20101221_AAAEEA roy_v_Page_73thm.jpg
5763d13b2f4d6bddc3b7ae16f03e1e83
8d5ca13a52d3b73c84ec582fefcbc7c8ece8695a
22057 F20101221_AAAEDL roy_v_Page_65.QC.jpg
cb84b98e9ea45ba46a073de6484877b9
bcf52a62b1c81d13bd63a47fc14717d162845d25
857711c19948c9dd030aa15d8b97e96f
a92a1cd21209f05a1bd837f5201b96d22438c4c7
28378 F20101221_AAAECY roy_v_Page_58.QC.jpg
8895e8edf7a50b0d294f7b939a6ca6be
2cac2739cb88176e0b17263d97ff5e04089e23dc
afd2b146bc1b7436292e485c147c7eef
8e33943e8d54133c900fb2fa4833ae798907a244
23227 F20101221_AAAEEB roy_v_Page_74.QC.jpg
f5e6c81a4dd24e674d34a43c02bc2e80
4975460be156b4536a37b9a09c0a3accf8e46830
20986 F20101221_AAAEDM roy_v_Page_66.QC.jpg
d34bfcca21b32cd18d1eac6499477fcc
43062bd51585b4d79da95885ee83274ab878c457
a2f315d0507267957d12586250bd6e31
2caee55e613bf83cde221f50a055a36996c7f4d5
6743 F20101221_AAAECZ roy_v_Page_58thm.jpg
647751ae469be95cc57f658c325a96fe
e94f4dcba6d9ab31e2f05ac2c54d0837d493bae9
4f7032cbc13a522657133e23397bbb0d
5d5e76944c6b410ee1b582fdbc7807690ffa2014
5993 F20101221_AAAEEC roy_v_Page_74thm.jpg
7f102e452bc9c9a76dd5ba93a51ec883
dcce7ba3c6a659e7311cb2084052914d0630b6ff
5756 F20101221_AAAEDN roy_v_Page_66thm.jpg
5445002616a7ec33aaaca5bd9f80860d
ff9e1fa422c80b854ed85bcdf23cc3bcde1d8fca
217bd308b01f71cba2e944d918696e7c
bb77d32fa7edf95fca9571a6c4031f4e7e6c6176
fb5292c5726e0fc1bc49bc3c7c29130f
7328bc57243bc1678a58c39bc9c9c0925e838dc5
18784 F20101221_AAAEED roy_v_Page_75.QC.jpg
742daea90a0cd9cf12bc8d285ed8bd59
6cde61f08acfbd8a9597734ed293cfe03cf7b58d
29197 F20101221_AAAEDO roy_v_Page_67.QC.jpg
cdf9268153dcaaef6af3de9963f5f5e2
d4f42dbab7f3edce8ac28012faa20b3a
96e7ed7eab642916b99dd47603995486c2e422a1
7521 F20101221_AAAEDP roy_v_Page_67thm.jpg
0eb4e091ca858f0de254a02968a7fdca
7b5279317b40491351b1dd208e0a25b5e6869373
4463801bdce80e7bccc4fc14337b676f
4ea8bab54c333be450c359f9a14b64dcbc4f2168
4dc549d9e339f28e6a412dc239dc1ee8
c7e2a3bb6835fd4207cdeac2ba5ef19504cb2d67
7316 F20101221_AAAEEE roy_v_Page_76thm.jpg
c9aa505c434a2711a700035f1a85192b
cf56d6fea8c3a98771acd04cef68474e0ea0a811
21434 F20101221_AAAEDQ roy_v_Page_68.QC.jpg
fd8e9f5a07ae37693d8b1807443af472
8ca16e045d227177f51acb908d462122d61a0f5f
ffb2ff8c579e553c2d5a33e34531bb76
f32305e4a60f3dd373afaa80186f347b
2d150f9f8c28e451778d57afec59701019bf9353
21561 F20101221_AAAEEF roy_v_Page_77.QC.jpg
8fa65b60da7032aae30284ed4326bb97347611fe
6124 F20101221_AAAEDR roy_v_Page_68thm.jpg
5d2c733fde49719695a844c69878e6ed4127d886
49445be5d6d7cda367dd716b6ee23950
b1bc41ab938de4f827e34f2237af1411b7ed1b92
5521 F20101221_AAAEEG roy_v_Page_77thm.jpg
e74d05a2bcbd710e9fec64ccb7e34936
6217a9645b83fb781d4c233bd035ec2740d56b3e
24975 F20101221_AAAEDS roy_v_Page_69.QC.jpg
15a3216390f4918be3ba746f49aa662c
1db0a034b3440082ac12df9a6f7680b49b8be7f8
e9a44a71834769e8d70170c8848c7f04
67a9f1b6cf8161db496b073a4c2073946a906dac
b4bde950266b5bcf3d895fcd1931f2a4
036cd103e70673f7050611d0d59fe9e2268c6d64
25043 F20101221_AAAEEH roy_v_Page_78.QC.jpg
2fd92c87a505c6cb2634c9a59a53eb85
a090165f8846d3669b90c179b38b47df77188f56
99273e0d38e9e868debbd7d0ab4c3e32
43aa4f668e3d9ceb975278b935592f8365568242
6173 F20101221_AAAEDT roy_v_Page_69thm.jpg
d11579a8a6fa776356048243452e1c0ef3a2ccc1
bd9efd5a8241defe24a934f284607bf3
0d2078dac0a0b24b739725092e99bec39d378f19
774dc7d91d0ae76959a771cbe14b4336
e7be99c703580baaba5194ee155f1a503ee824ff
6566 F20101221_AAAEEI roy_v_Page_78thm.jpg
2c0b9a07fd1d7ee146600ec4dcbd1682
c4eab3c6724aeccc0a6d242871327114df2cfba7
e342a9157914af7b17792e81bce2e154
292d0b23ec8a19702a987a0fe23e85bb311b5dc7
3958 F20101221_AAAEDU roy_v_Page_70thm.jpg
d7b43c84a9c8d36ff70227d6ce05e3c0
85660beb65029e02f043a640951a05d323fa5827
25188 F20101221_AAAEEJ roy_v_Page_79.QC.jpg
d64641deeb70a4c376e607670fd10c4a
bb904cefe4bb71ce3c85e23290afd2a35dbf667f
4d8546cbee92ec035a46a4510617b3ee9495e97e
27535 F20101221_AAAEDV roy_v_Page_71.QC.jpg
e98119abe5a36be66ee18b75ac8d0d0e
715c3674b551c74930160ec0eb841bbbc26f55dc
944691a9c381bc4bb43f9cd40e4175eb
053678eaa72c87ee1fca62d58aee173be4c52603
6480 F20101221_AAAEEK roy_v_Page_79thm.jpg
c6d2ccda6fb02a9ca75a3ec66e4a9d88
a267177fb0c4cb417f1deeb04ea313c5
4d4bf85ef2f7d5076613b65987382941d935e885
6932 F20101221_AAAEDW roy_v_Page_71thm.jpg
a6070e349492a621824a4f7f2a57cfda
e7380488ec04e21bc790887961534bfc19823d20
2c0804ec2cf5f61a14d9773cc31b6691
1779368a6043a261d46d5e30beb9f4f7f6a17793
23446 F20101221_AAAEFA roy_v_Page_89.QC.jpg
d35a59bce34daccf3b02944907d85ebe
252114f4b8d796c01d6fb23342d70e607ef70343
25696 F20101221_AAAEEL roy_v_Page_80.QC.jpg
f23f66179813966228f88e97fcfe157b
08a0a80054a9520fc6c3eaf768604a3304748aee
16367eebd21e3b63bc4163421e8eca98
17f9346c25663742f4e16a16715b25433dc58b9b
20153 F20101221_AAAEDX roy_v_Page_72.QC.jpg
ddf147f1b32fe43e5679368c436764b9
b685b0a9c553f5caf61aa9ed9d2c7dbd61dfe8ce
a69e9bd1b55a4599354aaaff41627f0e
32512 F20101221_AAAEFB roy_v_Page_90.QC.jpg
f93ace798b2bedc95618fdd80ae4a132
84e1169aa6779eef329e607a69d815543a9b1a30
28090 F20101221_AAAEEM roy_v_Page_81.QC.jpg
72a9fc14605d32aedb41c27978d8080c
9bce2ea53ac1e3e480143d57aa3c102e5a86c836
d4b44a31bcb563309f73871e4d3e6c43
35d1ef215d29c715307f898a9933c68fa023f348
5537 F20101221_AAAEDY roy_v_Page_72thm.jpg
b5079079157741a7a00c898d2dc5e83c
7821 F20101221_AAAEFC roy_v_Page_90thm.jpg
ffc10c26cbd871c488bdf9d01fed7ae5
ba6f0c4708904c258e22a678b9bd4e059e004dbb
6788 F20101221_AAAEEN roy_v_Page_81thm.jpg
288eb84c50f0c087b1dc88bea64b526e
d588c3af8e534dd1da8ba04951e947d4a1a913a5
bd9c08306d990c734ee2d00c869b78be82fb3742
28868 F20101221_AAAEDZ roy_v_Page_73.QC.jpg
e85664a4f8118fef9972150dc4a7ddea
6e2698d46e31402808d74be652f0d64322ff9f60
da096f430885ecb9660a66baf9e34e57
32680 F20101221_AAAEFD roy_v_Page_91.QC.jpg
417b7c954e1f18b64b87a631c0ee7c01ea9572a6
33136 F20101221_AAAEEO roy_v_Page_82.QC.jpg
8ca3c4abb6756125b5134b4bfb8e55e1
2932477992b3da278c79b126f559a9b1fddd047a
64f3d5e28bafc04de5f5a482eb2a086c
0e8b3a9139a21149c1c5b52b0672258ed8745654
33251 F20101221_AAAEFE roy_v_Page_92.QC.jpg
b569d355665d5b334ecb7c3c4ea857e3
462291e7be5a77b32a5080d7051bdaf986a50716
7740 F20101221_AAAEEP roy_v_Page_82thm.jpg
53a5e80e065e72e4dfa8527f7cbb1dba
75cebfa1a37c5b648ee85dbb49840cda549d8f88
9a180f8ea6bdbd9035167327c11aff0d
2bb0f6aef855824b8bb611cb39cd4d4f485fde0e
7bcb5929627ceea5d709e1e43fd2e214c538265b
23826 F20101221_AAAEEQ roy_v_Page_83.QC.jpg
8b7e977023bdfce6ccc0750c0b8367c0
8996b6c29607be81758eff3063f55095eae2813b
0b88b6222443a0dfa2a81f30da964abb
25aa8d90c42d95036c0f775cbee3fa7945942135
77b065a7130d1b271639b50ed65f6915
28dde5263d6eea12c7bb6596f848e3ee8b9981d0
8039 F20101221_AAAEFF roy_v_Page_92thm.jpg
24334f1c8b16b86079f71776727d278f
150a51226e60d3dd556cc6d8b4cfe167ee4b0111
7154 F20101221_AAAEER roy_v_Page_83thm.jpg
59ddf8063f45ee9326407568abdfe230
c15a482cc97a4fe5580b962e8e4eeaaa13a94e21
4e089bd3035c873cfff1fcc0a71af794
56f00deb016c63df357021b080244e465932279f
6cbaf2724cedf44804f621ef2c740953
69703fc867145e62ea64451f436f97eb53eecbdf
13624 F20101221_AAAEFG roy_v_Page_93.QC.jpg
8b483cee874bab5d806f177db749c0a65171112d
27463 F20101221_AAAEES roy_v_Page_84.QC.jpg
773d8f26abc99fc8e6524cd07b950b29
0cf997d0093f8f1a739151666e68ea1c5ac12fcc
27d2f1a906c03207af4a81de5329a987
a0f2732cd4c676aaa3249548530fd2b7fd9b3c4c
3468 F20101221_AAAEFH roy_v_Page_93thm.jpg
f54ea616ee6808db064344bd3b4760ea
f15120f10792a27731b198f6fe2f5b03e82fa426
6691 F20101221_AAAEET roy_v_Page_84thm.jpg
a50fd5a338dee86a630ef64eb408465c
ca3a77690fd621967446cf0471a7dd1fa13884f2
4dd1ae5bdf920c6ed45b00f0c9cf36ed
2809 F20101221_AAAEFI roy_v_Page_94thm.jpg
e65aba80025044b818c9f227a594a981
30072 F20101221_AAAEEU roy_v_Page_85.QC.jpg
3eb659fc00cbaeb92b4098d7c748378d
e904da0bd92acce442bab43c53542cea397281ea
09f1b316829187d3aa96b9c59ce26fd0
108073 F20101221_AAAEFJ UFE0022377_00001.mets
ff93ba8767a79635f953e3ee95b7c79c
c3c573cdc4a2d79cf966ff2aff7dee8ccb5c36b7
7345 F20101221_AAAEEV roy_v_Page_85thm.jpg
23c7dc6078035a2bc13087364d71bb0c
53d4843bb9502058e5662e5695f5d1c466f844f7
4a9966db16b7e4e1716cb705ceb110cd
685bb27d70780b34aab5eb45172a98dc3f688828
22345 F20101221_AAAEEW roy_v_Page_86.QC.jpg
098f207f86f0df19dffee20269e823bc
30045941a82f6430ccf4af19e4816c2b8bbd0124
560a5226740e5beaed2d4d7913354860
d8f2d1b067948c633d55b152228dd9ed7170e866
6129 F20101221_AAAEEX roy_v_Page_86thm.jpg
12b9bd1c426fa6edd576502ee1f2dfce4b650d4e
8815caaa1dd984f9564ef030b70ece69
7f0aec99ff4feb0ce863a6c7aa9faba36f427e41
5961 F20101221_AAAEEY roy_v_Page_87thm.jpg
b3d60899f41ca13ca4f57fa1b5596dda
5430d41c63d22051849d000f7ec666a1f081a2b0
455abbe4c9b8b48393c635e7e1a9f85528c60dc5
4028 F20101221_AAAEEZ roy_v_Page_88thm.jpg
1192a0530c676beab0e5e6e99b502cc0
986359eaef3d74d5d12fca21a382feba
2c07f2d50510b5a096c23b48956ac37db11a7644
f7bbacb08d14f5d5062b2d216a96fd89
ef422df3deb321a3a0678cd88c131aa4d7974e5e
929b4034e103bdafd39737d50e145a3a
91f1e84113684d43da22c62d0d07de14df457a05
e1951448a7a503085fd039a8b18d6a6c
4dc623ce72f0bce3b86b0c76b42a9cabf8c5fff1
f0e5ec3cb1cfe4f1844d903015b14ede
2ac5053b548d370faa3f8bfa942e604660347a5c
ecdcd993bcb58735eb4af4f28c6d7a79
a6bcea07a00e766ef13de5eb3f6647b3531f4ac4
56ddf78b8f5b17a00b72952fc8f575f4
1718e964a18d4cc4743a52e80870860edebda55b
d2cab3a0049ce6417d3ba017baf9dd2a819cf65a
08c1e45c4e0bccc9303bb903a8a828db
fa4726d7473ee381ea22a330a1c747976133f9fb
dab75a1ea49a897afae783f3d65dd6d7
63ea3178f9ce6e754ac618322e3eae60212c9569
9c7905cec83ab936f624f4967214454d
1f544a089afef595d5c94c5ee1b6059fe7662f37
06e58a07328ae10983c20110f1a273c3
1877570d9f303da1c0c6165dd448510b5714516b
0c2d03f5a0b6e6db03f639640831d13a
9ca793ffb4414f900c48566bde9599d1
44e1129d92c409dacf2302cffcdb9808f091b8ea
8d94d5038887885e095e0590677b62ae
6b712caf0d1791c41b1eabe0f1b1f8249ed754db
ab0e63a4be9b380e1ea6935a475bd85a
87ce2bc3964a3dccc663b40acce7be9995027ec1
6850150c60810f463a1ff5290c7148d8
406c70a0b2ffb15454436f36be14e1a5f6ec742b
60bf2be058404928c0f9b594e9624880
3b1af72d63b7fe8b4633cc7721ddeca4
159438c3412b246a443f6a27f36396bf8f98ced6
375dc1d53d7ff3ef9c1158275a6494e2
285b3753d3c6f9a12148081dbb00cbceec530019
ab3281cfee01764597a0eba001cd743d
70a950b32fa08478a433e769bb8b307e53876eb4
cf3f4fd4c624c86674abd2cfc5a1c0a5
22311d9b0011a220c8c79b76ef8a440010239277
5eef71058e0e824daaa5924118e1dd34
3159402a1119b7aa07736e91af41fa6f57c5b6b6
332d773b41424d1bcfc3f47c3e6df8f8
69035333eb8af68dfec7252667981b35ce6b256b
1b041254f35dc1de81115215fc8f659e
3cb2d24b46bec3ea831f8d1ccbf6a75e00d99391
4b026268dc44490f1c85b09c8314f389
9bdb2c0e581d94bdd73252afbdd30d24ded80c70
da8d03c99e88d6653541fd9455c87048
020ce34b57b13a6a1d2ed07d64c4f386
e9a0da6252d87f7e810362428982d89104185dd6
414665cc0f81f0401846fd91c9e9fdc8
5040a1b213d83864aa07834889041b3b41087504
295e46ccf6a83fb6be4e66f70d4e69b7
2a3ffaf65954f0eaa58d45466ac379683aebf363
4c9a2bb70ab1bfd7e9963ed733fe544e
b51555c7ab4fc4cddffe1884708f0046f0e97e09
b60a47a94c7dd5949be4701b18b4924e
5c56183ae86005e6c4ae5f40542365cf092fa7e9
7033d284b8cfb5fd3443fd500d3d240c
fec67d0ba58a8571d7eb786f0399852181679e83
af796cbe9dbed4afec5fd5ab708de0d9
ff69dcfbda74f7aed3037e22833a7c8e1abd69d6
4b6c344a5fa1dce1b46f8dc110648db1
6846fa205929f7a760a76101d9835bbbd3f3cde7
cdfbd7421ef686fa4c3691b71ce937cd
5c4d3df82c195a931baf9461b59cde2de61b8ea8
4ebe210c80f87e93080eee6fb0a4c819
849ff5509583518d881b96a9595774746c88723c
864980a1cc9c16054c5fc18488b3bd17
2f46c507a765199e8303ee3690d8fc0602edd04f
f594328b711281d33aab440733d89ff7
0e4293a37de8f1861fd81fc94cc30ec53b9dae4a
a16fa01c2c33855a9aefa13c9d67b719e035ec6a
cc99e1f891608aea3036e87b5a37ccd5
43c5d49779af108f42594911ca4e65af514dd388
8230f7f4aff6769602667a61cc7282ee
12ea6ba1e679d1a638b7385940e5ecbaf562fb0f
b496e099d34893be6a8184d7ae925552
3c93812c67615cc5520b37100c0d0c3c016b4af3
7b4c6a0069ea59ddb586db614f362cf6
2c5cf0a801c78d21631c78d2e84b99674e8ef212
e77d1c4fa050af6266ec21325b6a8b11
6c7feca2a94b6c5f3dd132379e13e05258e2a2f2
376662feb411d6a8035d9bef376fa02c
1d7794d97f279cc0123f1d177f755b4f184bfa1e
e2f3a2dc15c0a9ed1943b2d4209bb700
0fd63cb7dc948020758489decfc8f5a7e8c920be
9ded33a854ab23251b6c6200e229d3de
9d7d3bf61f0eab1fff02f42dc3bf0818cb47c0d3
eee145e554d6e2fa69f2338fd59bd5b5
840f624c245affda95f4628bed0ccf00b09a683d
7a8fdb5c72bb8cd35dea87b5e05fc1f3
e5c4306005995b65b6714e984ceda00bf04204fe
6575ea79d918417315695af784323f82
4641bf25efcd7247fc00613bcfa91b08e9e93d4a
0cbf6aaa3435e01ac8c2ec7dfe7c50a6
d2c07e4fb528edaf07f424c72a6b8dedfd84db27
a884d8805f62ddfa1aa602a318d2ed5b
889e58145bef61eba3ddec3d07d1365e
21af0c222cb518c0ed740c988c11c78e5ed73187
906e0e75a13c619be5e6bddb2ee99528
9739ca8d3eedcae2a859767ecf43ff5990742d3e
d22966a024e5dee685fb9b50eeab05e8
c7586b61013d5b0c68e8879c6bf5070f864846af
1e0566cdd4f88f0c0cfdfeee3dc0fc43
b245a5a176b0f4e89d0f06dc0eb0d58029a913f2
0a15dfbec6d02568c30afe2c28dcaf88
453f4f5f8bab4e612919fdeaeb0a22eb
42846781f2856da5ae922fe0c1d1cd8a9959956c
b15c77a9e68b7d7ab7e13947cdbe9cb7
fa388d9b2246590f2425f334f2d654cc
a01e360c3082d4820be6d87233ff6bb7ff123060
658bfc693151d41e8a5fd16059c9125b
6afdb395dd0612c6d3e2d8013928be758969ba68

THEORETICAL AND METHODOLOGICAL DEVELOPMENTS FOR MARK(OV
CHAIN MONTE CARLO ALGORITHMS FOR BAYESIAN REGRESSION

By
VIVEK(ANANDA ROY

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

2008

S2008 Vivekananda Roy

To my parents

ACKNOWLEDGMENTS

I extend my sincerest thanks to my advisor .Jint Hohert for his guidance throughout

my graduate study at University of Florida. I feel fortunate to have .Jint as my PhD

advisor. His guidance, help, enthusiasm, were all crucial to making this thesis take its

current shape. I'ni deeply grateful to hint for many other things, not the least for his

inspiring words in my hours of need.

I would also like to thank Professors Ben Bolker, Hani Doss and Brett Presnell for

agreeing to serve on my coninittee. I am particularly grateful to Professors Hani Doss and

Brett Presnell for being so kind to me over the past five years. I learned not only statistics

but also a lot about Emacs, LaTex and R front them.

I would also like to thank Professors Bob Dorazio, Malay Ghosh and Andrew

Rosalsky for sparing a lot of their valuable time on academic discussions with me and

giving me advice on several issues. I thank all my teachers from school, college and Indian

Statistical Institute whose dedication to teaching and quest for knowledge have inspired

me to pursue higher study.

Special thanks go to All .Ilv:adi and Parag for their friendship, care and support. I

have learned a lot about life in past five years front both of them. I owe deep gratitude to

.Jethinia whose care and affection I will never forget. I am thankful to Shuva whose love

and enthusiasm for niathentatics have ahr-l- .- inspired me.

I am indebted to many other people, mostly front my village, who guided me and

encouraged me during the formative years of my life: Arunda, Arunkaku, Bapida,

Bonikeshjethu, Budhujethu, Shashankajethu and my uncle.

Finally, I would like to thank my parents for ahr-l-~ .14eing a driving force in my life.

I often feel that whatever I have achieved is only due to nly parents' sacrifice, hard work

and honesty.

page

ACK(NOWLEDGMENTS .......... . .. .. 4

LIST OF TABLES ......... ... . 6

ABSTRACT ............ .......... .. 7

CHAPTER

1 INTRODUCTION ......... ... .. 9

2 MARKOGV CHAIN BACKGROUND . ..... .. 18

3 BAYESIAN PROBIT REGRESSION . ..... .. 24

3.1 Introduction .. .. ... . .. ..... .. 24
3.2 Geometric Convergence and CLTs for the AC Algorithm .. .. .. 27
3.3 Comparing the AC and PX-DA Algorithms .... .. 35
3.4 Consistent Estimators of Asymptotic Variances via Regeneration .. .. 37

4 BAYESIAN MULTIVARIATE REGRESSION .. .. .. .. 48

4.1 Introduction ......... . .. .. 48
4.2 Proof of Posterior Propriety ........ ... .. 50
4.3 The Algorithms ......... . .. 58
4.3.1 Data Augmentation ....... .. .. 58
4.3.2 Haar PX-DA Algorithm . ..... .. .. 60
4.4 Geometric Ergodicity of the Algorithms ... .. .. .. 61

5 SPECTRAL THEOREM AND ORDERING OF MARKOGV CHAINS .. .. 71

5.1 Spectral Theory for Normal Operators .... ... 71
5.2 Application of Spectral Theory to Markov C'!s I!us< ... .. .. 81

APPENDIX: CHEN AND SHAO'S CONDITIONS ... .. .. 89

REFERENCES ......._._.. ........_._.. 90

BIOGRAPHICAL SK(ETCH ....._._. .. .. 94

LIST OF TABLES

Table

page

3-1 Results based on R = 100 regfenerations ...... .. . 47

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

THEORETICAL AND METHODOLOGICAL DEVELOPMENTS FOR MARK(OV
CHAIN MONTE CARLO ALGORITHMS FOR BAYESIAN REGRESSION

By

Vivekananda Roy

August 2008

C'I I!r: James P. Hohert
Major: Statistics

I develop theoretical and methodological results for Markov chain Monte Carlo

(jl\C'l\C) algorithms for two different B li-. -1 Ia regression models. First, I consider a

profit regression problem in which }*\,...,Y*, are independent Bernoulli random variables

such that Pr(}* = 1) = #(:rf#) where :rs is a p-dimensional vector of known covariates

associated with }*, /3 is a p-dimensional vector of unknown regression coefficients and

#(-) denotes the standard normal distribution function. I study two frequently used

MC1| C algorithms for exploring the intractable posterior density that results when the

profit regression likelihood is combined with a flat prior on 73. These algorithms are

Albert and Chih's data augmentation algorithm and Liu and Wu's PX-DA algorithm.

I prove that both of these algorithms converge at a geometric rate, which ensures the

existence of central limit theorems (CLTs) for ergodic averages under a second moment

condition. While these two algorithms are essentially equivalent in terms of computational

complexity, I show that the PX-DA algorithm is theoretically more efficient in the sense

that the .I-i-ini!!lle'~ variance in the CLT under the PX-DA algorithm is no larger than

that under Albert and Chih's algorithm. A simple, consistent estimator of the .I-i-inia..'l~e

variance in the CLT is constructed using regeneration. As an illustration, I apply my

results to van Dyk and Meng's lupus data. In this particular example, the estimated

.I-i-inidul'lic relative efficiency of the PX-DA algorithm with respect to Albert and Chih's

algorithm is about 65, which demonstrates that huge gains in efficiency are possible by

using PX-DA.

Second, I consider multivariate regression models where the distribution of the errors

is a scale mixture of normals. Let xr denote the posterior density that results when the

likelihood of n observations from the corresponding regression model is combined with

the standard non-informative prior. I provide necessary and sufficient condition for the

propriety of the posterior distribution, xr. I develop two 1\C1| C algorithms that can he

used to explore the intractable density xr. These algorithms are the data augmentation

algorithm and the Haar PX-DA algorithm. I compare the two algorithms in terms of

efficiency ordering. I establish drift and minorization conditions to study the convergence

rates of these algorithms.

CHAPTER 1
INTRODUCTION

Realistic statistical modeling often leads to a complex, high-dimensional model that

precludes analytical, closed-form calculation which is required for statistical inference and

prediction. If we combine the complex model with a prior distribution on the unknown

parameters, as is done in B li-o -1 .Is Statistical analysis, the result is typically an intractable

posterior distribution of the model parameters given the observations. Suppose xr(0|y)

is the posterior density of the p x 1 vector of unknown model parameters, 8, given

the observations, y. In B li-, -i Ia inference, we are often interested in evaluating the

expectation of some function, my- f, with respect to the posterior density xr, i.e., we want

to know

Ex f= f 0)x(|y) O .(1-1)

Because the density xr(0|y) is a complicated function, closed-form calculation of the

above integral is generally impossible. We assume that the above integral exists and is

finite. Since Ex f can not be evaluated analytically, we use either deterministic numerical

integration techniques or simulation based methods to get an approximate value of (1-1).

Before delving into these computational methods, we provide two motivating examples.

In both of these examples, statistical modelling results in an intractable posterior density

making explicit closed-form calculation of the corresponding posterior expectations

impossible.

In problems involving toxicity tests and '.1. .-- li experiments, the responses are

often binary since what is observed is whether the subject is dead or whether a tumor

has appeared. A popular method of analyzing binary data is through B li-o -1 .Is analysis

with a probit link function. Suppose that we observe n independent Bernoulli random
variables, Yi,...,Y,, and we assume that Pr( = 1) = #(xf ) where x4 is ap xn 1 vectonr of

known covariates associated with ~, is a p x 1 vector of unknown regression coefficients

and #(-) denotes the standard normal distribution function. For ye {0(, 1}", that is,

Y = (Y1,... ,ys~) and y, E {0,1}), we have

i= 1

If we use a flat prior on p, the marginal density of the data takes the form

It is not obvious whether cy (y) < 00. We address this issue in chapter 3. Assuming

city) < 00, the posterior density of p takes the following form

i= 1

Clearly, the posterior density xr(P | y) is too complicated to allow explicit closed-form

calculation of posterior expectations of functions of P.

It has long been known that heavy-tailed error distributions are often required

when modelling financial data. Specific scale mixtures of normal distributions can be

used for modelling heavy-tailed data. Our second example is a B li- Ion multivariate

regression model where the distribution of the errors is a scale mixture of normals.

Suppose Y1, Y2, *, n are d-dimensional random vectors (e.g. returns on some assets)

satisfying the linear regression model

S= PXi + Ei

where p is the k x d matrix of unknown regression coefficients, xi's are k x 1 vectors of

known explanatory variables and we assume that, conditional on the positive definite

matrix E, the d-variate error vectorS E1, ,E are independently and identically

distributed with common density

fH6d E)= exp E ETZ-1E dUH(5
o(2xr) |ECa

where H(-) is the distribution function of a non-negative random variable. The density,

fH, clearly, iS a multivariate scale mixture of normals. The density fH can be made

heavy-tailed by choosing H appropriately.

We can rewrite the above regression model as

Y = Xp +E

where Y = (Yi, .. ., Y,)T is the a x d matrix of observations, X = (xl, xa2, *, )T iS the

a x k matrix of covariates and E = (E, .. En) iS the a x d matrix of error variables. The

likelihood function for this regression model is given by

f y|, ) (x)|E exp 6 -0x 2 H6

If we consider the standard noninformative prior on (P, E), i.e., if we assume that the prior
d~l
density, xr(P, E) oc |E| 2 the posterior density takes the following form

C2 1 x "

C2(y i i=1 0~ (2xr) |E|-I 2

where c2 9) is the marginal density of y given by

d(d+1)
where W CR 2W is the set of d x d positive definite matrices. In chapter 4, we

provide necessary and sufficient conditions for c2 9) < OO. As in the previous example,

posterior expectations with respect to the posterior density, xr(P, E | y), are not available in

closed-form.

We now discuss different computational methods that can be used to approximate

(1-1). These computational methods are broadly of two types, namely, numerical

integration methods and simulation based methods. If the dimension, p, is not large,

numerical integration techniques can be efficiently used to obtain a good approximation of

(1-1). But, as p increases, numerical integration techniques become less and less efficient

because of the well known problem called the curse of I.:Is:.: ,:r.......lU;, In this dissertation,

we consider simulation based methods to estimate the posterior expectations.

An alternative to numerical integration is to estimate (1-1) by Monte Carlo sampling.

Monte Carlo integration requires drawing iid samples X*,X*,...,XA_, from xr(-) and

then using the sample mean
m-1
fm f(X J ),
j=0
to estimate the expectation (population mean) in (1-1). The justification of Monte Carlo

methods comes from the strong law of large numbers (SLLN~), which guanrantees that fm

converges almost surely to Ex f as m tends to infinity. So, Ex f can be well approximated

by fm provided the sample size, m71, is large enough. We often know r only up to its

normalizing constant, i.e., usually, ,, xr(0|y) de is unknown. In that case, we can use

rejection sampling methods Robert and Casella [38, C'!s Ilter 2.3] to obtain an iid sample

from xr. Rather than giving details about different Monte Carlo methods, we now address

an important issue that the experimenters ahr-l- .- face -i.e., how to choose the agi

sample size, m?

How large a sample size is sufficient is a subjective matter. It depends on how much

error we are willing to accept in the approximation. One way to measure the accuracy

in the approximation is by the width of a 95' confidence interval for Ex f. A confidence

interval for Ex f can be obtained using the central limit theorem (CLT) for the estimator

fm. If f has a finite second moment with respect to r, i.e., if Exf2 < 00, then by the
classical central limit theorem we have

fl7, Ex f 1 N(0,v2) aS n 4 OO

where v2 = Exf2 (Exrf)2. An .l-i-mptotic 95' confidence interval for Ex f is given by

Im + 2sm,/~ where s,, is a~ strongly con~sistent, estim~a~tor of v2" giVen? by
m-1

i=0

Since the sample size, m, is under our control, the main benefit of calculating the

confidence interval is to determine whether the Monte Carlo sample size we choose is

large enough. In practice, one draws a random sample of size M~ from xr for some finite

n~umber M! anld constructs th~e confidence interval fMI + 28M/2il. If the Ilength? of the

resulting confidence interval, 4sM/ lz, Seems to be satisfactory, then one stops the

sampling and reports fM RS an eStimate of Ex f. On the contrary, if the .I-i-mptotic 95' .

confidence interval is deemed too wide, then M~ can be increased appropriately and further

simulation can be carried out until the desired level of accuracy is achieved. Of course,

in the latter case, if we know beforehand the precision, e, that we want to achieve then

we can use SM/ RS a pilot estimate of v to calculate an approximate sample size, namely,

(28M @l~2, that we need.

In practice, making iid draws from xr might not be feasible. For example, in the probit

regression model that we mentioned, it is difficult to make iid draws from the posterior

density, xr(Ply), especially when the dimension, p, is large. Similarly, for the multivariate

regression model that we discussed, it is problematic to produce a useful Monte Carlo

method to simulate from the posterior density, xr(P, Ely).

Surprisingly, it is straightforward to construct a Markov chain with stationary

distribution xr even when direct simulation from xr is impossible. As explained in the next

paragraph, it turns out that it is indeed possible to approximate (1-1) by simulating a

Markov chain with stationary distribution xr. This is the basic principle of Markov chain

Monte Carlo (ifi'lC) method. The most general algorithm for producing Markov chains

with arbitrary stationary distribution xr is the Metropolis-Hastings (il-H) algorithm. A

simple introduction to the M-H algorithm is given in Chib and Greenberg [7]. Another

widely used MCijlC algorithm is the Gibbs Sampler [4]. Suppose the p-dimensional vector

8 in (1-1) can be written as 8 = (01, 82, ). The simplest Gibbs sampler (but, not

the general Gibbs sampler) requires one to be able to simulate from all univariate full

conditional densities of xr i.e., it is required to simulate from the conditional distributions,

s| {0;, j / i} for i = 1, 2,..., p. It is also possible to create a hybrid algorithm which uses

different versions of M-H algorithm together with Gibbs sampler to construct a Markov

chain with stationary distribution xr. As our discussion so-----~ -r- there is a plethora

of Markov chains with stationary distribution xr. In order to choose between MCijlC

algorithms, we need an ordering of Markov chains having the same stationary distribution

xr. In C'!s Ilters 2 and 5, we describe different partial orderings of Markov chains.

Let {Xjy}Ro denote the Markov chain associated with an MCijLC algorithm that is

used to explore xr. If {Xjy}R is Harris ergodic (defined in C'!s Ilter 2), the ergodic theorem

implies that, no matter what the distribution of the starting value, Xo,
m-1
fm := ~f(X,)
j=0

is a strongly consistent estimator of Ex f, i.e., fm Ex f almost surely as m c o. So,

the ergodic theorem (like the SLLN in iid case) ensures that Ex f can be approximated

by running a well behaved Markov chain for sufficiently large number of iterations. In

practice, as in iid case, one simulates the chain for a finite number of iterations, ;?i M~',

and reports fM;r aS the estimate of Ex f. Suppose there is an associated central limit

theorem (CLT) given by

fm x f (0, .2) S 0 00 (1-)

and that we have a consistent estimator of o.2, Call it ~2. So, we can compute an

.I-i-inidlli'lc standard error for fM;r, Which is given by &/ ~l. As in the iid case, the

.I-i-inidlli'lc 95'. confidence interval given by fMr & 29-/21 can then be used to decide

whether there is any need for further simulation.

As -II_0-r-- -1. in the previous paragraph, establishing a central limit theorem for a

Markov chain is essential in order to put MC10L on equal footing with iid sampling.

Unfortunately, unlike in classical Monte Carlo methods, the finite second moment

condition i.e., E, f2 < 00 does not insure a CLT for f,~. In addition, the Harris ergodicity

which establishes the strong consistency of foz is not enough to guarantee that (1-2)

holds. It generally requires rigorous analysis of the Markov chain {Xj;}Ro in order to

prove that CLT holds for f,,. There are several v- .--s Of establishing the CLT in (1-2).

These approaches can he broadly divided into two categories. One approach is based on

probabilistic (convergence rate) aI, l1i--;-; of the Markov chain. We give a brief description

of these techniques in C'!s Ilter 2. The other approach exploits results from functional

analysis (see Chapter 5).

Another difficulty in constructing the confidence interval, fn/r & 2o-/ W1 is that even

when there is a CLT, finding a simple, consistent estimator of the .I-i-mptotic variance,

0.2, can he challenging due to the dependence among the random variables in the Markov

chain. Mykland, Tierney, and Yu [:33] show that when CLT exists, regenerative simulation

(R S) methods can he used to construct a consistent estimator of o.2 by uncovering the

regenerative properties of the Markov chain (Section :3.4). The regenerative simulation

technique basically breaks the whole Markov chain up into iid pieces (tours) by keeping

track of the regeneration times. Then, standard iid theory can he used to analyze the

.I-i-ini!!lle'~ behavior of the ergodic average, f,,2 and thus a simple, consistent estimator

of .l-i-!Injdull- variance is obtained. It might not he easy to implement the RS method in

practice. There are other methods like batch mean and spectral methods which are easier

to employ to estimate the .I-i-inidull-; 1 variance (Jones et al. [21] and the references cited

therein). The advantage of using RS method is that it is on stronger theoretical footing

than the other methods.

We now provide a brief overview of the four remaining chapters of this dissertation.

In the next chapter, we review some results from general state space Markov chain

theory. In particular, we mention sufficient conditions for Markov chain CLT and provide

a partial ordering of Markov chains based on their performance in the central limit

theorem.

In (I Ilpter 3, we study two MC \!C algorithms that are frequently used for exploring

the posterior density, xr(Ply), that we mentioned before in the context of the probit

regression example. These algorithms are Albert and Chib's [1993] data augmentation

algorithm and Liu and Wu's [1999] PX-DA algorithm. We study the convergence rate

of these algorithms and prove the existence of central limit theorems (CLTs) for ergodic

averages under a second moment condition. We compare these two algorithms and show

that the PX-DA algorithm should ah-- .--s be used since it is more efficient than the

other algorithm in the sense of having smaller .I-i-mptotic variance in the central limit

theorem (CLT). A simple, consistent estimator of the .I-i-inidllicl~ variance in the CLT is

constructed using regenerative simulation methods.

In C'!s Ilter 4, we consider B li-, -i Ia multivariate regression models where the

distribution of the errors is a scale mixture of normals. We noticed before that if

the standard noninformative prior is used on the parameters (P, E), then posterior

expectations with respect to the corresponding posterior density, xr(P, Ely), are not

available in closed-form. We develop two MC1| C algorithms that can be used to explore

the density xr(P, Ely). These algorithms are the data augmentation algorithm and the

Haar PX-DA algorithm. We compare the two algorithms and study their converge rates.

We also provide necessary and sufficient conditions for the propriety of the posterior

density, xr(P, Ely).

While in ChI Ilpters 3 and 4, we used probabilistic techniques to 2.1, lli-. .. different

MC1| C algorithms, it is possible to take a functional analytic approach to study and

compare different Markov chains. In ChI Ilpter 5, we give a brief overview of some results

from functional 2.!! l1i--;- In particular, we discuss the spectral theorem for bounded,

normal operators on Hilbert space. We show how these results of functional analysis can

be used to study Markov chains.

CHAPTER 2
MARK(OV CHAIN BACKGROUND

Let A = {Ami~o denote a time-homogeneous discrete-time Markov chain on a

general state space X equipped with a countably generated o--algebra B(X). Let Pm(x, A)
be the m-step Markov transition function associated with A for m = 1, 2, 3,. .. So,

Pm(x, A) denotes the probability that the Markov chain at x will be in the set A after m

steps (transitions), that is, for x E X, Ae B (X) and le {0(, 1, 2, ...},

Pm(x, A) = Pr(Am+i E A | At = x)

When m = 1, we simply denote the one step Markov transition function by P(x, A) and

for m = 2, 3, .. ., Pm(x, A) is defined iteratively by

Pm~x A) P~, dyPm-(y, A).

A probability measure xr on B(X) is called an invariant probability measure for A, if, for all

measurable sets A,

xr(A) = P,, AxIdx)

Note that,

CP2(x, A)a(dx:) = SP(x.)P, d)P(y, )=~x = (dy)P(y,l A)= (A).

Similarly, we can show that f Pm(x, A) x (dx) = x (A) for m = 1, 2, 3,. .. So, if Ao ~ xr,

then Am ~ xr for all m and A is stath-.: o a,~r in distribution.

Let L2(;,) be the vetor space of real-valued, measurable functions on X that are

square-integrable with respect to xr, i.e.,

The Markov chain A is said to be reversible with respect to xr if for all functions f, g e

L2(r

fC S (yl~gr)P(x, dy)r(dx) = f~S I(x)g(y)P(x, dy)i(dxl).

If we take g(x) = 1 in the above equation, we get

f S~(y)P(x, dy~~h)=XS(dx) =f(x)Ix dy)i(dx) = f l~(x:)xd)

i.e., xr is invariant for A. So, if a Markov chain is reversible with respect to xr then xr is
invariant for the chain.

Suppose that p is a non-trivial, o--finite measure on X. The Markov chain A is
called p-irreducible if for each x E X and each set A with p(A) > 0, there exists an

m EN := {1, 2,. .. }, which may depend on x and A, such that Pm(x, A) > 0. In words,

the chain is p-irreducible if every set with positive p measure is accessible from every

point in the state space. The measure p is called an irr..I;. 09ii;, measure for A. As in
M.~ in and Tweedie [30, Section 4.2], when we ;?i "A is ~-irreducible" we mean that A\

is p--irreducible for some p and that is a maximal irre 7 09ii;, measure for A. Two

properties of maximal irreducibility measures that will be used in the sequel are (i) if p is
an irreducibility measure and is a maximal irreducibility measure, then p is absolutely

continuous with respect to (denoted > p), and (ii) a maximal irreducibility measure is

unique up to equivalence, i.e., if I1 and 2a are both maximal irreducibility measures, then

#1 > #2 and 2a > 1 (denoted 21 a

The p--irreducible Markov chain A is ap~eriodic if there do not exist an integer d > 2

and disjoint subsets Ao, Al,. ., Ad- _c X with p~(Ao) > 0, such that for all i = 0, 1, .. d- 1

and all x E Ai,

P(x, Aj) = 1 for j = i + 1(mod d).

Suppose A is ~-irreducible and define B+(X) = {A E B(X) : ~(A) > 0}. The Markov

chain A is called Harris recurrent if for all Ae B +(X),

Pr (Am EA i.o. | Ao = x) = 1 for all x E X.

The Markov chain A is called Harris ergodic if it is ~-irreducible, periodic and Harris

recurrent. The Harris ergodicity of a Markov chain is often easy to verify in practice and it

implies that, for every x E X,

||Pm(x, -) jT(-) | 0a m co ., ,

where ||Pm(x, -) xr(-)|| denotes the total variation distance between the probability

measures Pm"(x, -) and 7i(-), i.e., the supremumn over measurable A of Pm"(x, A) x(A) .

However, the Harris ergodicity tells us nothing about the rate at which this convergence

takes place. If it takes place at a geometric rate, then A is said to be geome/0.:. ellol ergodic.

More precisely, the Harris ergodic Markov chain A is geometrically ergodic if there exists a

constant p E [0, 1) and a function M : X [0, 00) such that for any x E X and any me N ,

|| P~x,-) -x () | < Mx) m .(2-1)

We now describe methods that are used to prove the geometric ergodicity of a Markov

chain. One method of proving that A is geometrically ergodic is by establishing drift

and minorization conditions. There are several v-wsi~ of doing this (11. i-n and Tweedie

[31], Rosenthal [45], Roberts and Tweedie [40]). Here, we describe a method based on

Rosenthal's [1995] work.

A drift condition holds if for some function V : X [0, 00),

PV < AV + L

for some A E [0, 1) and some L < 00, where (PV)(x) = f, V(y)P(x, dy). The function V is
often called a drift function.

An associated minorization condition holds if for some probability measure Q(-) on

B(X) and some E > 0 we have

P(x, A) > EQ(A) VX E C and VA E B(X)

where C := {x E X : V(x) < l} with I being any number larger than 2L/(1 A).

Rosenthal's [1995] Theorem 12 shows that the above drift and minorization conditions,

together, imply that A is geometrically ergodic. In Chapter 4, we employ drift and

minorization conditions to prove the geometric ergodicity of the data augmentation

algorithm used in B li-, Io multivariate Student's t regression problem.

One advantage of proving geometric ergodicity of A by establishing the above drift

and minorization conditions is that using Rosenthal's [1995] Theorem 12, we also can

calculate an upper bound of M~(x)pm in (2-1). This upper bound can be used to compute

an appropriate burn-in period ( Jones and Hobert [23], Marchev and Hobert [28]). There

are other methods of proving geometric ergodicity of a Markov chain that do not provide

any quantitative bound of M~(x)pm in (2-1). We describe one such method now.

We will assume that X is equipped with a locally compact, separable, metrizable

topology with B(X) as the Borel o--field. A function V : X [0, 00) is said to be

unbounded of compact sets if for every y > 0, the level set {x : V(x) < y} is compact.

The Markov chain A is said to be a Feller chain if, for any open set Oe B (X), P(-, O) is

a lower-semicontinuous function. The following proposition is a special case of M.~ i-n and

Tweedle's [1993] Lemma 15.2.8.

Proposition 1. Sup~pose that the Harris ergodic M~arkov chain A is a Feller chain.

Supplose further that the support of a maximal irrech;.-.l..7. ;, measure has non-tii pl ;

interior. If for some V :X [0, 00) that is unbounded of compact sets

PV < AV + L

for some A E [0, 1) and some L < 00, then the M~arkov chain, A, is geomen,... ellol ergodic.

In C'!s Ilter 3, we apply Proposition 1 to establish geometric ergodicity of MC \!C

algorithms used in B li-, -i Ia probit regression problem. Hobert and G. o;r [15] emploi-. I

Proposition 1 to establish the geometric ergodicity of Gibbs samplers associated with
B li-, -i Ia hierarchical random effects models.

Notice that, unlike Proposition 1, the drift condition in Rosenthal's [1995] Theorem

12 does not require the drift function, V, to be unbounded off compact sets. Also,

Rosenthal's [1995] Theorem 12 does not need A to be a Feller chain.

The driving force behind MCijlC is the ergodic theorem, which is simply a version

of the strong law that holds for well-behaved Markov chains, e.g., Harris ergodic Markov

chains. Indeed, suppose that f : X R I is such that |x Ifldx < 00 and define Exf =

f dx The theergoic teore .v that the average fm = m-] CE f(As) converges

almost surely to Ex f no matter what the distribution of Ao. This justifies our use of fm as

an estimator of Ex f. We will ;?i that there is a CLT for fm if there exists a a2 E (0, 00)

such that, as m 00o,

L(f m Ex f) iN~(0, a"2

As explained in OsI Ilpter 1, CLTs are the basis for .I-i-inidllicl~ standard errors, which can

be used to ascertain how large a sample is required to estimate Ex f. Unfortunately, while

the Harris ergodicity of a Markov chain does imply that the ergodic theorem holds, this

is not enough to guarantee the existence of CLTs. However, if A is geometrically ergodic

and reversible with respect to xr, then the CLT holds for every f such that S, f2d~ o

that is, for every f E L2(;T) [41]. (For more on the CLT in MCil C, see C'I I>. and Geyer

[5], Mira and Geyer [32], Jones [20] and Jones et al. [21].) For a thorough development of

general state space Markov chain theory, see Nummelin [34] and M.~ i-n and Tweedie [30].
Robert and" Rosentha^l [43] provides a concise, self-contained description on general state

space Markov chains (also see Tierney [49]).

As mentioned in OsI Ilpter 1, for a given distribution function, xr, there are large

number of MCil C algorithms with stationary distribution xr. One way to order these

algorithms is based on their performance in CLT. Note that the .I-i-md.l!lle~; variance, O.2

in (1-2) depends both on the function f and the particular MC1| C algorithm that we

are using. Suppose P and Q he the Markov transition functions corresponding to two

different MCijlC algorithms with stationary distribution xr. Let us denote O.2 for these

two algorithms by n( f, P) and v( f, Q) respectively. Assume, both v( f, P) and v( f, Q) are

finite. Then if we are interested in calculating E, f, we prefer the Markov chain P over Q

if v(f, P) < v(f, Q) provided the two chains are equivalent in terms of simulation effort.

On the other hand, if we do not assume any prior knowledge about the function whose

expectation we want to evaluate, we need a uniform ordering as below.

Definition 1. [SE] If P and Q
Markov chains with intericent Igo~~~l.:.:sl..;i measure x., then P is better than Q in the

e~ff~.,.. I ordering written P FE Q. if c(f, P) < 'U(f, Q) for every f E L2(,)

In ChI Ilpter 3 and ChI Ilpter 4, we order different MCijlC algorithms in terms of

efficiency ordering.

CHAPTER 3
BAYESIAN PROBIT REGRESSION

3.1 Introduction

Suppose that Yi,...,Y, are independent Bernoulli random variables such that

Pr(~ = 1) = #(xTP) where xi is a p x 1 vector of known covariates associated with ~, P

is a p x 1 vector of unknown regression coefficients and #(-) denotes the standard normal

distribution function. For ye {(, 1}", that is, y = (yl,..., y,) and yi E {0,1}), we have

i= 1

A popular method of making inferences about P is through a B li-o -1 .I analysis with a flat

prior on p. Define the marginal density of the data as

C'I. in and Shao [6] provide necessary and sufficient conditions on y and {xi}", for

city) < 00 and these conditions are stated explicitly in the Appendix. When these

conditions hold, the posterior density of P is well defined (i.e., proper) and is given by

Unfortunately, the posterior density xr(P | y) is intractable in the sense that expectations

with respect to it, which are required for B li-, -i Ia inference, cannot be computed in closed

form. Moreover, as we mentioned in (I Ilpter 1, classical Monte Carlo methods based

on independent and identically distributed (iid) samples are difficult to apply when the

dimension, p, is large. These difficulties spurred the development of Markov chain Monte

Carlo methods for exploring xr(Ply). The first of these was Albert and Chib's [1993] data

augmentation algorithm, which we now describe.
Let X denote the, \x p design matrix whose ith row is x { and, for z = (zi, ..., z,)T E

RW", let p = P(z) = (XTX)-1XTz. Also, let TN(p, is2, w) denote a normal distribution with

mean p and variance is2 that is truncated to be positive if w = 1 and negative if w = 0.

Albert and Chib's algorithm (henceforth, the "AC als..i s~I lIn~ ) simulates a Markov chain

whose invariant density is xr(P | y). A single iteration uses the current state P to produce

the new state p' through the following two steps:

(i) Draw zi, .., z, independently with ze ~ TN(x'P, 1, yi)

(ii) Draw p' ~ N,( i(z), (XTX)-1)
Albert and Chib [1] has been referenced over 350 times, which shows that the AC

algorithm and its variants have been widely applied and studied.

The PX-DA algorithm of Liu and Wu [27] is a modified version of the AC algorithm

that also simulates a Markov chain whose invariant density is xr(P | y). A single iteration of

the PX-DA algorithm entails the following three steps:

(i) Draw zi, .., z, independently with ze ~ TN(x'P, 1, yi)

(ii)~ ~ ~ ~ 1 Dra g2l ~ am liz x(XI1XI)I; 1Xz)2 and set z' =(z,.,g,
(iii) Draw ii' ~ Np( j(zj), (X X)-1)
Note that the first and third steps of the PX-DA algorithm are the same as the two steps

of the AC algorithm so, no matter what the dimension of P, the difference between the

AC and PX-DA algorithms is just a single draw from the univariate gamma distribution.

For typical values of n and p, the effort required to make this extra univariate draw

is insignificant relative to the total amount of computation needed to perform one

iteration of the AC algorithm. Thus, the two algorithms are basically equivalent from

a computational standpoint. However, Liu and Wu [27] and van Dyk and Meng [51] both

provide considerable empirical evidence that autocorrelations die down much faster under

PX-DA than under AC, which so-----~ -r- that the PX-DA algorithm "mixes f I-I. I than the

AC algorithm. (Liu and Wu [27] also established a theoretical result along these lines see

the proof of our Corollary 1.)

Suppose we require the posterior expectation of f (P) given y, i.e., we want to know

assuming this integral exists and is finite. Let {@}Rj"o denote the Markov chain associated

with either the AC or PX-DA algorithm. We later show in this chapter that {@}Rj"o is

Harris ergfodic. So the ergfodic theorem implies that, no matter what the distribution of

the starting value, Po,
m-1

j=0
is a strongly consistent estimator of E [ f(n) | y] ; that is, fm, E [f ( ) | y] almost surelyv

as m 00o. As defined in C'!s Ilter 2, we ;?i- that there is a CLT for fm if there exists a

0.2 E (0, 00) such that, as m 00o,

ImT, E [ f() | y]) N(0, 02) aS ,n ix OO .-1)

As explained in OsI Ilpter 1, establishing the central limit theorem for fm is crucial to make

honest statistical inference based on {py }Ro. We know that one way to ensure CLT in

(3-1) is by establishing geometric ergodicity of {@}Ro. In this chapter, we prove that the
Markov chains underlying the AC and PX-DA algorithms both converge at a geometric

rate which implies that the CLT in (3 1) holds for every fe L2(7i(lj Iy)); that is, for

every f such that fy, f2(P)xr(PIy)dp < 00. We also establish that PX-DA is theoretically
more efficient than AC in the sense that the .I-i-mptotic variance in the CLT under the

PX-DA algorithm is no larger than that under the AC algorithm. Regenerative methods

are used to construct a simple, consistent estimator of the .I-i-ing d o)tic variance in the CLT.

As an illustration, we apply our results to van Dyk and Meng's [2001] lupus data. In this

particular example, the estimated .I-i-inidllicl~ relative efficiency of the PX-DA algorithm

with respect to the AC algorithm is about 65. Hence, even though the AC and PX-DA

algorithms are essentially equivalent in terms of computational complexity, huge gains in

efficiency are possible by using PX-DA.

The remainder of this chapter is organized as follows. Results that the AC and

PX-DA algorithms are geometrically ergodic appear in Sections 3.2 and 3.3, respectively.

In Section 3.4 we derive results that allow for the consistent estimation of .li~!!llh d ic

variances via regenerative simulation.

3.2 Geometric Convergence and CLTs for the AC Algorithm

We begin with a brief derivation of the AC algorithm. Let RW = (0, 00), R_

(-oo, 0], z = (zi, .., z,)T E R" and let #(v; p, x2) denote the N(p, x2) density function

evaluated at the point ve R Consider the function from RW"x RW" R W given by

x(P,v z )=Is z) )ys a zeIo v)(zi; xf f )

where, as usual, IA (. iS the indicator function of the set A. Note that

i= 1

and hence, xr(P, z | y) can be viewed as a joint density in (p, z) whose marginal is the

target density xr(P | y). This joint density is usually motivated as follows. Let Z1, ..., Z,

be independent random variables with Zi ~ N(xTP, 1). If we define = Ipg (Zi), then

Yi, .. ,Y, are independent Bernoulli random variables with Pr(~ = 1) = # (xT f). The Zi's

can therefore be thought of as latent variables (or missing data) and xr(P, z | y) represents

the posterior density of (p, z) given y under a flat prior on p. The AC algorithm is simply

a data augmentation algorithm (or two-variable Gibbs sampler) based on the joint density

xr(P, z | y). Indeed, a straightforward calculation reveals that

| z y ~ N ( (X'X)X )-1

and conditional on (P y), Z1, .. ,Z, are independent with

Zi | 79, y ~ TN(.r {/, 1, yi).

If we denote the current state of the Markov chain as /9 and the next state as /9', then

the Markov transition density of the AC algorithm is given by

Note that k(P | 79) xr(P I| y) = k(P I | 7') xr(P | y) for all /3, /' E RIW'; i.e., k(P | 79) is reversible

with respect to r(fi | y). It follows ininediately that the posterior density is the invariant

density for the Markov chain, or, in symbols,

Let K(-, ) denote the Markov transition function corresponding to the AC algorithm;

that is, for /9 E RIW and a measurable set A,4

The corresponding m-step Markov transition function is denoted by K'"(ft,4). We now

show that the Markov chain driven by k(fi' | 79) is Harris ergodic.

Let p- denote Lebesgue measure on RIW. Several nice properties follow front the fact

that K(fi, -) has a (strictly positive) density with respect to p. Indeed, if 7t(,) > 0, then

K((P, A) > 0 for all /9 E RIW'; i.e., it is possible to get front any point 79 E RIW to the set ,4 in

one step. This implies that the AC algorithm is ys-irreducible and aperiodic.

In order to establish Harris recurrence, we must introduce the notion of harmonic

functions. A function b : RIW" R is called harmonic for K if h(fi) = (Kh)(fi) for all

/9 E RI>. One method of establishing Harris recurrence is to show that every bounded

harmonic function is constant [34, Theorem 3.8]. Suppose b is a bounded, harmonic

function. Since the AC algorithm is (-irreducible and has an invariant probability

distribution xr(P | y), it is recurrent, which in turn implies that & is constant ~-a.e. [34,

Proposition 3.13]. Thus, there exists a set NV with p(NV) = 0 such that h( ) = c for all

pe N Now, for any pe R P~, we have

which implies that h c. It follows that the AC algorithm is Harris recurrent.

We have now shown that the Markov chain corresponding to the AC algorithm is

Harris ergodic and thus from C'!s Ilter 2 it follows that ergodic theorem holds for it. The

following theorem is the main result of this section.

Theorem 1. The M~arkov chain on RW with transition 1:.: l 0 k (P' | P) (that is, the

Markov chain underlying the AC rlly.>rithm) is geomen,... ellol ergodic.

Proof. We will show that the AC algorithm satisfies the hypothesis of Proposition 1.

We have shown that AC algorithm is p-irreducible and periodic, where p denote the

Lebesgue measure on RP". So, if is a maximal irreducibility measure for the Markov

chain underlying the AC algorithm, then > p. Conversely, if p(A) = 0, then Km(P, A)=

0 for all pe R P? and all me N which implies that ~(A) = 0 and it follows that p > ~.

Hence, pa -. Since the support of p obviously has non-empty interior, it follows that the

support of a maximal irreducibility measure for the AC algorithm has non-empty interior.

We now demonstrate that the Markov chain associated with the AC algorithm is a

Feller chain. Let P and O denote a point and an open set in R ", respectively. Assume

that {#1}"7,, is a deterministicc) sequence in RP~ with p& / such that pt i as

m c o. Two applications of Fatou's Lemma in conjunction with the fact that xr~z|4, y) is

continuous in yield

lim inf K (PA,O) > lim inf k (p'| m ~) dp

= lim inf x( R* | lnz ) z
J to mm

> x~fi' z, y) int inf ~rx~a|7,aU

=K(P O) ,

and hence K(-, ) is a lower-senticontinuous function. Hence, the Markov chain corresponding

to the AC algorithm is a Feller chain.

We apply Proposition 1 with drift function V(fi) = (X/S)T(X/S). Recall that X is
assumed to have full colunin rank, p, and hence XTX is positive definite. Thus, for each

,* > 0, the set

P liE RI'? : V <) I 19=( E RIW : pTXTXp IS <

is compact so the function V is unbounded off compact sets. Now, note that

(KV) (p> ~ i)kf'7)p

= EE [V(:i') r-, ] i?-Y y

where, as the notation elo----- -r- the expectations in the last two lines are with respect

to the conditional densities xr(f' | x, y) and xr~x | 7, y). Recall that xr(f' | x, y) is a

p-dintensional normal density and xr~x | 7, y) is a product of truncated nornials. Evaluating
the inside expectation, we have

E [V(P ') zy] = E [(P ') X X/S' z,y]

=tr(X X(X X)-1) +: XX(X X)-1(X X) (X X)-1X :

=p + X X(X X)-1X :

< p+X z,

where tr(-) denotes trace of a matrix and the inequality follows from the fact that

z I-X X(XX -1TX z>0

for all z E R"n. We now have that

EE ;[V(I') z, ] ;), y < Ep +z z y =p+ E [z: | ,y].
i= 1

Standard results for the truncated normal distribution [19] imply that if U ~ TN((, 1, 1)

then,

E(U2) ~ 2

where #(-) with only a single argument denotes the standard normal density function; that

is, #(v) is equivalent to ~(v; 0, 1). Similarly, if U ~ TN((, 1, 0) then,

E(U2) ~ 2

It follows that

1 +( x pj ) ( if y = 0 .
A more compact way of expressing this is as follows:

SLiE [zf2 | ,y]= 1+(re-i) 4)2 ,/,Tn (3-2)

where I, is defined in the Appendix. Hence, we have

(KVT)(i) = E E[V;(P') zly] i,y) i= i

Recall that the goal is to show that (KV)(P) < AV(P) + L for all PR E W. It follows from

(3-3) that (KV)(0) < p + n. We now concentrate on P eRI \ {0}.
We begin by constructing a partition of the set RW \ {0} using the a hyperplanes

defined by wTP = 0. For a positive integer m, define Nm = {1, 2,..., m}. Let

Al, A2,..., Aan denote all the subsets of No, and, for each j E IT_ define a corresponding

subset of p-dimensional Euclidean space as follows:

Sj = {p3E R"\{0} :wfgT~<0fo rall i EAjand wfg7P >0fo raill i EAj}

where Aj denotes the complement of Aj, that is, Aj = N, \ Aj. Note that

the Sj are disjoint,

U zlSj = Rw \ {0}, and
some of the Sj may be empty.

We now show that if Sj is nonempty, then so are Aj and Aj. Suppose that Sj / 0 and

fix p E Sj. Since the conditions of Proposition 5 are in force, there exist strictly positive

constants {ai }", such that

Therefore ,

aimT + a~lTP + ---+ I, ,tiP = 0. (3-4)

The matrix X has full column rank p, and hence 0 < pTXTXp = CE (xTP)2

E = /ET 72. Thus, there exists an ie N such that wTP / 0 and, since all the ai are
strictly positive, (3-4) implies that there must also exist an i' / i such that I, 4 n r

have opposite signs. Thus, Aj and Aj are both nonempty. Now define C = { je E : T

0}. For each j E C, define
/ieA T n2 ieA.~ zT n2

and

xj = sup Ry (P) E [0, 1].

In the following calculation, we will utilize a couple of facts concerning the so-called

Mill's ratio. First, when n > 0, a (u)/(1 #(u)) > u2 [11, p.175]. Also, it is clear that if

we define

M=sup
ue(-oo,o] 1 (u) '(U
then M~E(0, 00).

Fix j E C. It follows from (3-3) and the results concerning Mills ratio that for all

SE Sj, we have

i=1

< p+ n + (,, )2

1- P)(wP)
ieA ~(~i

iEAj

(KV) (0)

-` (w 4p)2
iEj

(I,.TP)~(WTP)
1-~(WTP)

ie:

iEAj

< p+ n+ (w, P)2 + nM
i=1

=p + ~n(M + 1) + (w 4Tij)2

= xvp +nM+ 1)+R()(L)

where L := p + n(M~ + 1). Therefore, since s~cS,

(K V) (P) < AV (P)+ L ,

where

A := max Xj
jec

Hence, it suffices to show that As < 1 for all j E C.

Again, fix j EC and note that for le R ,, Rj(10)

Rj(P) which means that Rj(P)

depends on P only through p's direction and not on its distance from the origin. Thus,

xj = sup Ry (P = sup Ry (P < sup Ry (P) ,
P6S3 PES3 PES:*"

RIW \ {0}, it follows that

where

Sf* = {p E RP : ||4|| = 1 and wffT < 0 for all i E Aj and wff'i > 0 for all ie Ay } ,

and

Sf* =I {p n:|4|=1adwf o l lA n wf>0frali s

Now sice j*is a compact set in RP" and Rj(P) is a continuous function on Sj**, we know

that

sup R ( ) = Ry ( ) for some p E Sj**

Assume that p E Sj* is sulch that R,( ) =1, that is,

-2 -2

This implies that CE ,l(ll'n O) Again, there exist strictly positive constants

al, a2, to Such that

a~wim + a~wp + --+ a,wT =0o.
But we already, know, that wf#T = 0 for ll i e Aj, and hence it must be the case that

Howecver, w{~i < 0 for a~ll i eA, as a Sf*. This combined with th~e fact thait as a~re all
strictly, positiven shows that w{,T = for all i E Aj. Hence, we have identified a nonzero

such that

ITr = 0 for all ie N,

But this contradicts the fact that W has full column rank. Therefore, we have established

that

sup Ry (P) <1 ,

which implies that As < 1. Therefore, A < 1 and the proof is complete.

Together with the results of Roberts and Rosenthal [41], Theorem 1 implies that

the AC algorithm has a CLT for every f e L2(i( l g. Ill Order to use this theory to

calculate standard errors, we require a consistent estimator of the .-i-mptotic variance, o.2

This topic will be addressed in Section 3.4. In the next section we show that geometric

ergodicity of the AC algorithm implies that of the PX-DA algorithm and that PX-DA is
at least as good as AC in terms of performance in the CLT.

3.3 Comparing the AC and PX-DA Algorithms

The Markov transition density of the PX-DA algorithm can be written as

k*(p'l | )=xp ',yRz z)>z|4 )d

where R(z, dz') is the Markov transition function induced by Step 2 of the algorithm

that takes z z = (gzi,..., gz,)T. It is straightforward to show that the Markov

chain driven by k* is Harris ergodic. Hobert and Marchev [17] provide results that

can be used to compare different data augmentation algorithms (in terms of efficiency

and convergence rate). In order to establish that their results are applicable in our

analysis of k*, we now show that R(z, dz') admits a certain "group representation." Let

"xz | y) = fy, xr(P, z | y) dp. A simple calculation reveals that

|XTX|-4 exp ( z (I H)z/2}
cy~ly (y 2x

where H = X(XTX)-1XT. Let G be the multiplicative group R,+ where group

composition is defined as multiplication; i.e., for gl, g2 E G, gl o g2 = 192. The identity
element is e = 1 and g-l /.Telf-armaueo su~g dg/g where dg

denot~es L~ebesgue: measure: onl R Let~ GJ act on the left of RW" through component-wise

multiplication; that is, if geG c andu z e R", then gz = (gzi,..., gz,). With the left group

action defined in this way, it is easy to see that Lebesgue measure on RW" is relatively left

invariant with multiplier X(9) = 9"; i.e.,

J~S R ("ll J Rild

for all geG c andu all int~egrable functions h : R's R I. (See ChI Ilpters 1 & 2 of Eaton [10]

for background on left group actions and multipliers.) Let Z denote the subset of RW" in

which x lives; i.e., Z is the Cartesian product of n half-lines (RW and RW_), where the ith

component is R,+ if yi = 1 and RW_ if yi = 0. Fix x E Z. It is easy to see that Step 2
of the PX-DA algorithm is equivalent to the transition x i gx where y is drawn from a
distribution:- on G having density function

X(9);g ( y)v4(d9) 9;z-ls(gz y)dyg (c'( H):) ZT- -g' (IHiz/2 d
.fe X(9)~rg ( y)v4(d9) .1 9;z-li(gx | y)dy 2(it-2)/2F(n/2) eg.

Furthermore, S, (g)xr(gx | y);4(dy) is positive for all x E Z and finite for almost all: E Z.

Consequently, we may now appeal to several of the results in Hohert and Marchev [17].

First, their Proposition :3 shows that R(x, dr') is reversible with respect to xr~x | y) and it

follows that k*(f3'|/3) is reversible with respect to xr(P3| y). We now use the fact that the

AC algorithm is geometrically ergodic to establish that the PX-DA algorithm enjoys this

property as well.

Corollary 1. The M~arkov chain on RP" with transition I7:. um / k* (f' | /3) (that ise. the

M~arkov chain underlying the PX-DA rlly~.:,thm) is geomtl,... ril;i ergodic.

Proof. Define

Let K and K'* denote the Mlarkiov operators on L (r(n3 y)) associated with the Markiov

chains underlying the AC and PX-DA algorithms, respectively [25, :32]. Denote the norms

of these operators by ||K|| and ||K*||. In general, a reversible, Harris ergodic Markov

chain is geometrically ergodic if and only if the norm of the associated Markov operator

is less than 1 [41, 44]. By Theorem 1, the AC algorithm is geometrically ergodic and

consequently ||K|| < 1. But LiU and Wu [27] show that ||K*|| < ||K|| [see also 17, Theorem

4] and hence ||K*|| < 1, which implies that the PX-DA algorithm is also geometrically

ergodic. O

We have now shown that the Markov chains underlying the AC and PX-DA

algorithms are both reversible and geometrically ergodic and hence both have CLTs

for all f e L2(r(p |)). We now use another result from Hobert and Marchev [1'7] to show

that the PX-DA algorithm is at least as efficient as the AC algorithm.

Theorem 2. Fix f e L2(P Iy).I a~,ka a~,k d6801 the UGnritRCeS in the CLT f~or

the AC and PX-DA rlly.., /thms,, i~ 1.: 0 .l;, then 0 < o)~,k* a~,k
Proof. The result follows immediately from Hobert and Marchev's [2008] Theorem 4. O

In order to put our theoretical results to use in practice to compute valid .I-noi- nd l'lc

standard errors, we require a consistent estimator of the .I-i-mptotic variance and this is

the subject of the next section.

3.4 Consistent Estimators of Asymptotic Variances via Regeneration

We begin with the AC algorithm. Instead of considering the Markov chain on RW

driven by k(P' | p), we consider the joint chain on RW"x RW" with Markov transition density

given by

The Markov chain dlefined by k, which we denote by {pyj, zy}Ro,1 has invariant dlensity

xr(P, z | y) and satisfies the usual regularity conditions. It may seem to the reader more

natural to ulse the Malrkov transition density k (z', p' | z, 4)= p'|zyxz| ,)

for the joint chain. We have discussed this issue in Remark 3. Of course, the marginal

chain {@}Rj"o has the Markov transition density k(P' | P) no matter which version of the

Markov transition density we choose for the joint chain. While this is obvious for the

chain corresponding k, it can be easily shown for k by considering two consecutive steps

of the joint chain. Let k(P'|4) be the Markov transition density of the {@},t"o chain
corresponding k. Suppose, two consecutive steps of the joint chain are (P, z) and (P', z').

Then ,

=. k(p'|z, P)k(z|4)dz

where the conditional densities in the third equality are obtained from the Markov

transition density, k, of the joint chain. The de-initializing arguments of Roberts and

Rosenthal [42] can be used to show that the joint chain {pj,zy},t"o inherits geometric

ergodicity from its marginal chain {pj },to Note that {#y, zy }=,o is the chain that is
actually simulated when the AC algorithm is run (we just ignore the zys).

Suppose we can find a function a : RW x RW" [0, 1], whose expectation with respect

to xr(P, z | y) is strictly positive, and a probability density d(P', z') on RW"x RW" such that for

all (P', z'), (p, z) ERI? x R"n, we have

k (P', z | 4 z) > S (r, z) d(P', z') (3-6)

This is called a minorization condition [22, 30, 43] and it can be used to introduce

regenerations into the Markov chain driven by k. These regenerations are the key to
constructing a simple, consistent estimator of the variance in the CLT. After explaining

exactly how this is done, we will identify s and d for both AC and PX-DA.

Equation (3-6) allows us to rewrite k as the following two-component mixture density

k (p', z' | 4, ) = s (p, z)d(P', z') + (1 a (#, z)) r (P' z' | 4, z) (3-7)

where r is the so-called residual density defined as

k (P',z' | 4, z) Stp,z)d(P', z')
1 s(P, z)

when s(P, z) < 1 (and defined arbitrarily when s(P, z) = 1). Instead of simulating

the Markov chain {pj,zy}Rj"o in the usual way that alternates between draws from

"xz | p, y) and xr(P | z, y), we could simulate the chain using the mixture representation

(3-7) as follows. Suppose the current state is (pj, zy) = (p, z). First, we draw 6j ~

Brcnoulli(s(P, z)). Then? if ri, = 1, wei draw (iy I:, zy) from~ dl, anld if by = 0, we draw~n

(pj 1, zy41) from the residual density. The (random) times at which by = 1 correspond

to regenerations in the sense that the process probabilistically restarts itself at the next

iteration. More specifically, suppose we start by drawing (Po, zo) ~ d. Then every

time by = 1, we have (pj 1, zy 1) ~ d so the process is, in effect, starting over again.

Furthermore, the "tours" taken by the chain in between these embedded regeneration

times are iid, which means that standard iid theory can be used to analyze the .I-i- pind icl~

behavior of ergodic averages, thereby circumventing the difficulties associated with

analyzing averages of dependent random variables. For more details and simple examples,

see Mykland et al. [33] and Hobert et al. [16].

In practice, we can even avoid having to draw from r (which can be problematic)

simply by doing things in a slightly different order. Indeed, given the current state

(pj, zy) = (P, z), we draw (pj 1, zy 1) in the usual way (that is, by drawing from x~z | p, y)

and xr(P | z, y)) after which we "fill in" a value for by by drawing from the conditional

distribution of 6j given (pj, z) and (pj 1, zy 1), which is just a Bernoulli distribution with

success probability given by

We now describe exactly how these supplemental Bernoulli draws are used to construct a

consistent, estima~tor of E [f(p) | y] a~s well as a con~sistentt estim~a~tor of the correspon~ding

.I-imph!lli(c variance.

Suppose the Markov chain is to be run for R regenerations (or tours); that is, we

start by drawing (Po, zo) ~ d and we stop the simulation the Rth time that a by = 1.

Let 0 = -ro < -rTi r~ 72 < R be the (random) regeneration times; that is,

-te = min {j > -re-1 : by _l = 1} for te { 1, 2, .. ., R}. The total length of the simulation, -rR,

is random. Let NI~, N2.., ** R be the (random) lengths of the tours; i.e., Nst = -rt T-rt1

and define

N\ote that the (Ns;, St) pairs are iid. The strongly c~onsisten testimnator of E [ f (P) | y] is

S1

j=0

whee 3= R 2=1 Stan N= 2<1 N. Because the Markov chain driven by k is

geometrically ergodic, the results in Hobert et al. [16] are applicable and imply that, as

long a~s there exists anl au > 0 such that E[| f(P)|2+ | y] < 00, then1

d fr-E~f0)|y N(072 9)

as R c o. (Note that the requirement of a finite 2 + a~ moment is a bit stronger than the

second moment condition discussed earlier.) The main benefit of using regeneration is the

existence of a simple, consistent estimator of y2, Which takes the form

2 i~ t=1St i

Remark 1. The CLT in (3-9) is lkib~l;, different from the CLT discussed earlier, which

takes the form

Hobert et al. [16] explain that the two CLTs are retlated by the equation 2 =" E[s(P,z,) | y]O.2

Remark 2. A further ral.. tr,:y.: of using regeneration to calculate standard errors is that

the starting distribution is prescribed to be d(P,zx) so that burn-in is a non-issue.

We now derive a minorization condition for the AC algorithm using the "distinguished

In .11.1 technique introduced in Mykland et al1. [33]. First, note t~hat k (p', z' | 4,z) does not,

depend on P and as a consequence, neither will our function s. Fix a distinguished point

z, E R"n and let D be an p-dimensional hyper-rectangle defined by D = D x x D,

where Di = [ci, di] and ci < di for all i = 1, 2, .. ,p. Now note that

=s(z) d( ', z')

where

s (z)= E inf and d (P', z') = 1 rx (' | p', y) x (P' | ze, y) ID( ') (3-10)

and

J~~,( RP""( J.vi (l~"~( RnJD
Clearly, d(P', z') is a probability density on RW x RW". All that is required to apply the

regenerative method described above is the ability to draw from the density d (to start

the simulation) and the ability to calculate fl in (3-8). Making a draw from d( ', z') can

be done sequentially by first drawing P' from the truncated density E-1~(P / X* Y) D P)

(which does not require the value of E) and then drawing z' from xr~z' | P', y).

-1X z, I

exp zX(XX-Xz

exp z X(XTX)-] X z

We now provide a closed form expression for s(z), which in turn will give a closed

form expression for the success probability r7 First,

x(4 | z, y) = 1 exp
(2x)9|XTX|-

(z)) .

21

li(z)) XTX (4

where P(z)

S(z)

(XTX)- XTz Thus,

x(P | z, y)
E inf

E inf'
06Dexp (z-,'~x(,) XT z) 0XTX (z,)]

exp'

z' X(X7X)

z,) Xp)

exp c gIi =]

(tii )

(ti) + d Iti p

exp

SzT X(XT X)

'X1X

i~onf exp(z

where tT = (z z,)TX. Therefore, the success probability rl in (3-8) becomes

=~j mf) IDp~l j+1)
s#xD d5 ~, xjy) j1j1 +1 9

= mf ID j+1i

~(~lpj~ex)~pj~ ct Is (@ ) y I

exp -zXXX-Xz
( )x Xp I

ex

- -~~~

j+1)

where ti)T 1 z -z)TX and py 1,4 is the ith element of the vector 4 ~. Note that, rl

depends only on pja1 and zy. Also, note that, E is not required to calculate rl.

Notice that there is a chance for regeneration only when the P component enters the

p-dimensional rectangle D. This so__~-1-;- making D large. However, increasing D too

much will lead to very small values of rl. Hence, there is a trade-off between the size of D

and the magnitude of the success probability, rl.

Modifying a computer program that runs the AC algorithm so that it simulates

the regenerative process is quite simple. Since code for simulating from xr( '|z, y) and

xr~z' | p', y) is already available, it is straightforward to write code to simulate from d. All

that remains is a small amount of code to calculate rl and compare it to a Uniform(0, 1)

after each iteration of the AC algorithm.

Remark 3. A minorization condition can also be obtained for k using the "distinguished

1p~...ni" technique. However, in this case, the z component has to enter an n-dimensional

i ),l

exp -T zf X(TX)~-X zy

exp [i I =) 1 I /(j) ) 1 6"ID j+

r LItri:ll~. before a regeneration is possible. In most applications, a is much Irlary ? than p,

and when n is Iary,l the <.-?.;?7.17~l;, that all n components of z -:::l,,tll~r, .... ;, le; enter their

assigned interval is ini'~..a ll;i so small that the rll' .-r:thm is of no practical use. Moreover,
as mentioned above, this problem cannot be solved .: cl~;, by I,,;;.;. .9l the intervals 'ary ,I .

Regeneration can also be used in conjunction with the PX-DA algorithm. Indeed, we
now show that a simple modification of our minorization condition for the AC algorithm

yields a minorization condition for the PX-DA algorithm. The Markov transition density
of the PX-DA algorithm can be rewritten as

k*(p' | 4)=xp z yhg|zxz| ,y gd

where h(g | z) is the density in (3-5). As before, instead of working directly with k*, we

consider the joint chain on RW"x RW" x RW with Markov transition density given by

k*(p', (z',g') P, (z, g)) = [h(g' | ')(z' | ', y)] x(P |gzy).

Let {#y, (zy, gj)} o denote the Markiov chain corresponding to k*. Since T(4 | y) is the
invariant density for k*(P'|4), we have

=( x)p | gz~~iy)L [hg| i~z|y) g z

and from this it follows that xr(P | z,: y) [h(g | z)xr(z | )] is the invariant density for

{pj, (zy, g))}Ro. As before, considering two consecutive steps of the joint chain, it
can be shown that the marginal chain {pjy}Ro has Markov transition density k*. It is
straightforward to show that the joint chain is Harris ergodic, and, as before, the chain

associated with k* inherits geometric ergodicity from its marginal chain {Pjy}Ro. We can

get a minorization condition for k* as follows

k*(p, (', g) |4, (, g) =[h.(g' |')7i(z' | /j', y)] (p' | gz, y)

cs( Is (t ) did Ig (t )

exp -~ zTX(XX-Xz

s(gjzy) d*(,By 1, (zy 1, gj 1))

.~ x(6gzyy (,Bl |zej, y)

expi gzX (XTX)-1X z p

exp z X (X X)-1X z,

exp g3 zX(X)Xz

=exp= ct *()Ipg (t () +dit*j Ip_ (t*j ) t ,y4ID;j1

j+1>

where t*C F) = (yzi z)7X. TheCoremI 2 states that the .l-i Intl'l' d ic va~rian~ce in? th~e CLTi

for the PX-DA algorithm is no larger than that for the AC algorithm, i.e, O < o-f,k*I

o-f~k < OO for all f e L2 (7(,6 |y)). However, we know from Remark 1 that the regenerative
method is based on a slightly different CLT whose .I-i-inidllicl~ variance has an extra factor

involving the small function, namely E(s() | y), from the minorization condition. Although

the small fumetions in the two minorization conditions that we derived for the AC and

PX-DA algorithms are slightly different, E(s(.) | y) remains exactly the same as shown

>l~ o( in x(,6 |)( gz, y) hg|z)xz|,,y) (6|zyID;

=s(gz)d*(,', (z',g'/)),

where the function s(-) is as defined before and

d* (,', (z', g')) = [h(g' | z')r(z' |,6', y)]x(,6' ze,y/)I7D ;6)

where E is alSo the same as before. Hence, for the PX-DA algorithm the function rl

becomes

i1

below

=~~g sIz~~ | y)z|y)1

i= 1
=~~~2 cos sz) '(I H)z') g"2- e-z' (I-H)z'/2

x exp, z' (I H)z'/2g2)C --n R {} i R {} i

1+ (z') (I H)7( z')l~21J (I~T(1~~X
2(n-2)/2F(n/2) (Z' g" ex (- '(I- )z/g2 /

where z-' = gz. Clonseque~ntly, if E[| f~(P)|2+ < OO FOT Some a? > 0 and if r ,k and

qj,k* denote the variances in the regenerative CLT for the AC: and PX-DA algorithms,

respectively, then 0 < yf~k* I ,k <00. Hence, PX-DA remains more efficient than AC in
the regenerative context.

We end this section with an illustration of our results using van Dyk and Meng's

[2001] lupus data, which consists of triples (yi, Xl, Xi2), i = 1,... ,55, where xsi and Zi2

are covariates indicating the levels of certain antibodies in the ith individual and yi is an

indicator for latent membranous lupus nephritis (1 for presence and 0 for absence). van

Dyk and Meng [51] considered the model

with a flat prior on p. We used a linear program (that is described in the Appendix) to

verify that C'I. in and Shao's [2000] necessary and sufficient conditions for propriety are

satisfied in this case.

In order to implement the regenerative method, we had to choose the distinguished

point z, as well as the sets [ci, di]. We ran the PX-DA algorithm for an initial 20,000

iterations starting from the maximum likelihood estimate of P given by P = (-1.778, 4.374, 2.428).

We took the distinguished point to be the average value of z over this initial run. For

is { 0, 1, 2}, let pi and as denote the sample mean and sample standard deviation of the

#is over this initial run. Wle set De [iA 0.09 si, A+ 0.09 s,]. (The factor 0.09 s, was

chosen by trial and error.)

While generating Gamma and multivariate normal random variables is straightforward,

we need an efficient algorithm for generating truncated normal random variables. We

used the accept-reject algorithm of Robert [39] to generate one-sided truncated normal

random variables. We ran AC and PX-DA for R = 100 regenerations each. This took

1,223,576 iterations for AC and 1,256,677 iterations for PX-DA. We used the simulations

to estimate the posterior expectations of the regression parameters and the results are

shown in Table 3-1. (Results in C'I. i. and Shao [6] imply that there exists a~ > 0 such that

E [|py |2+" ly < OO for je {0(, 1, 2}).) It is strikiing that the estimated .I-i-ing~l .i ic variances
for the AC algorithm are all at least 65 times as '609.-~ as the corresponding values for

the PX-DA algorithm. These estimates -II__- -r that, in this particular example, the AC

algorithm requires about 65 times as many iterations as the PX-DA algorithm to achieve

the same level of precision. (We actually repeated the entire experiment seven times and

th~e estimates of 7 ,k 7 ,k* ran~ged betweenl 40 anld 145.)

Table 3-1. Results based on R = 100 regenerations

AC Algorithm PX-DA Algorithm
Parameter estimnate s.e. (~fi;fk leStimnate s.e. (;jf;k* ,k
Po -3.060 0.097 -3.018 0.012 6;6.6;
P1 7.005 0.190 6.916 0.023 66.9
2a 4.037 0.121 3.982 0.015 63.1

CHAPTER 4
BAYESIAN MULTIVARIATE REGRESSION

4.1 Introduction

Suppose Yi, Y2, are d-dimensional random vectors satisfying the linear

regression model

where p is the k x d matrix of unknown regression coefficients, xi's are k x 1 vectors of

known explanatory variables and we assume that, conditional on the positive definite

matrix E, the d-variate error vectorS E1, ,E are independently and identically

distributed with common density

fH E) = exp E E-1E d1H(5) (4-2)

where H(-) is the distribution function of a non-negative random variable. The density

fH is a multivariate scale mixture of normals and it belongs to the class of elliptically

symmetric distributions. The density fH can be made heavy-tailed by choosing H
appropriately. For example, when H is a Gamma(" ,\ )distribu~tion fulncition, fH becomes

the multivariate Student's t density with v > 0 degrees of freedom i.e., the density of El

becomes proportional to [1 +-FT- v"E 2-6 II !!v of the results in this chapter a~re

specific to the multivariate Student's t regression model.

We can rewrite the regression model in (4-1) as

y = Xp +E

where y = (yl, .., y,)T is the a x d matrix of observations, X = (xl, x2., )T iS the

a x k matrix of covariates and E = (E, .. En) iS the a x d matrix of error variables. The

likelihood function for the regression model in (4-1) is given by

We consider the standard noninformative prior on (P, E), i.e., we assume that xr( E) oc

|E| 2 The posterior density takes the following form

where c2 9) is the marginal density of y given by

d(d+1)
where W CR 2W is the set of d x d positive definite matrices. Fernandez and Steel [12]

proved that c2 9) < OO if and only if a > d + k. In section 4.2, we give an alternative

proof of the posterior propriety. A byproduct of our proof is a method of exact sampling

from xr(P, Ely) in the particular case when n is exactly d + k. Throughout this chapter

we assume that n > d + k. We also assume that the covariate matrix, X, is of full

column rank i.e., r (X) = k. The posterior density in (4-3) is intractable in the sense that

posterior expectations are not available in closed-form. Also, our experience shows that

it is difficult to develop a useful procedure for making i.i.d. draws from xr(P, Ely). In this

chapter we focus on MC1| C methods for exploring the posterior density in (4-3).

We develop a data augmentation (DA) algorithm for xr( Ely) in section 4.3.1. It

has been noticed in the literature that the standard DA algorithm often suffers from slow

convergence [51]. Empirical and theoretical studies have shown that alternative algorithms

that are modified versions of the standard DA algorithm such as the Haar PX-DA

algorithm and the ;;;.,97, l.:1r augmentation algorithm often provide huge improvement over

the standard DA algorithm (Liu and Wu [27], van Dyk and Meng [51], Roy and Hobert

[46], Hobert and Marchev [17]). In section 4.3.2, we develop the Haar PX-DA algorithm

for the posterior density in (4-3).

We then specialize to the case when the errors, Ei S, have a Student's t distribution
i.e., the mixing distribution, H(-), in (4-2) is a Gamma(, ") disotribu~tion fu~nction.W

prove that in this case, under certain conditions, both the DA and the Haar PX-DA

algorithms converge at a geometric rate. As mentioned in Chapter 2, the geometric

ergodicity of a Markov chain guarantees the existence of central limit theorem. Using

results from Hobert and Marchev [17] we also conclude that the Haar PX-DA algorithm is

at least as efficient as the data augmentation algorithm in the sense that the .-i-mptotic

variances in the central theorem under Haar PX-DA algorithm are never larger than those

under the DA algorithm. Some of these results are generalizations of results from van

Dyk and Meng [51] and Marchev and Hobert [28] who considered the special case where

there are no covariates in the regression model (4-1), i.e., X = (1, 1,...,1)T and H is
Gamma(" \ ) distribu~tion fuinction

The rest of this chapter is laid out as follows. In the next section, we prove the

propriety of the posterior distribution. In section 4.3, we describe the DA and the Haar

PX-DA algorithms. In the last section, we compare the two algorithms and prove that

both the algorithms converge at a geometric rate.

4.2 Proof of Posterior Propriety

In this section we address the propriety of the posterior distribution. In particular, we

prove that xr(P, Ely) is a proper density for almost all y if and only if a > d + k.

Theorem 3. Let y be Lebesgue measure on R d". The posterior I1:.:; 0 Hr(P, |y) is

proper for y- almost all y if and only if a > d + k.

Proof. We want to show that

for y- almost all y if and only if a > d + k. Recall that the likelihood function f(y|4, E)

is itself an integral with respect to the distribution function H(-). We now use this fact

to rewrite c2 9) aS aHOther integral which is defined on a larger space and is easier to

handle. In order to simplify notation, we assume that H(-) has a pdf h(-) (with respect

to Lebesgue measure on R,+). We introduce latent data q = (ql, q2, > ) Such that,

conditional on (i, E), { (y;-, qj)) n are iid pairs satisfying the relations

and qj|#,C E h(-).

If we denote the joint conditional density of y and q hy fly, q|/3, E), then by (4-1) and

(4-2) it follows that

I

f(U 9 /S, E)49

f(UO91j,, / )f(4jP9) / E4j

f(U /SE). (4-5)

So we get

JW JR

JW JR L JR" J

R" WRf(y q|/S, E)x(/3, E)d/S dEdq

where the last equality follows by Fubini s theorem. Note that, the integrand in the right

hand side of the above equation is given by

II
f~ ~ ~~~n (y, q|/3 EYj i ) exp
(2xr) z|E|lz =1

2 )=1

I
/(UO/S E)

Let Q be an ax n diagonal matrix with diagonal elements (- -,i~ -). Then,

j=1

j=1

j= 1

j=1

=tr E-p 0 -p p pyQ

(X Q- X)-1XTQ-ly. The following definition of the

Arnold [2, C'!s Ilter 17]

Yj)TC-I(PTXj

where, = (XTQ-1X)-1 and p =

matrix normal distribution is from

Definition 2. [] Suppose a~ is an mxr matrix;, E and T are mxm and rxr non-negative

7. I;,../.: matrices. We .rt;, that Z has a matrix: normal distribution with parameters a~,

and T if Z is an m x r random matrix: having moment generating function given by

and we write Z ~ Nm,r (a, E, T) In this case we have E(Z) = a~. Moreover, if E and T are

positive 7. It...:I~ matrices then Z has the following I. ,.-.:1/ function

1 1
.fz(z) =ep -t{ ( ) ( ) .
(2xr)mr/21 r/2 rm/2 2

Since r(X) = k, it follows that X Q-1X and hence R is a p.d. matrix. Thus,

(2r)9 | |0| 2I2I exp -
(2xr) 2 |E|?

Itr E -y -p || hq)E -
2

d+(n-k)+ 12

(2xr) 2 2 I = ( i
We now assume that n > d + k and will show that c2 9) < Xo. We f Tst show that when

n > d + k, y Q- y <"0- ft is a p.d. matrix with probability one (with respect to the

Lebesgue measure on R~d"). Notice that

yTQ-ly rgS-1q = Y -1Y YT -1X(X Q-1X)-1X Q-1Y

= g~ yi~x Q- I-Q ( X)-X'iX Q-Q- (4-7)

and since I Q-2X(XTQ- X)- X Q-2 is an idenipotent matrix, it implies that

yTQ-ly ,TS2-if is a positive senli-definite matrix. Now we prove that y Q-ly ,rg2-11
is a p.d. matrix by showing that |yTQ-ly 79S-1<| / 0 (with probability one). Let A be

the n x (d + k) augmented matrix (X : y). Then,

A'Q-1 A = Q X

X Q-1X X Q-1Y

y Q-1X Y Q-1Y

Therefore ,

|A'Q- A| = |X Q- X||Y Q-1 y -YTQ- X(XTQ- X) XTQ- Y|

= X Q-1X||y Q-l 79S-11

= |C-ll 7 -l 79S-19|. (4-8)

Since r(X) = k, we know that X Q-1X is a p.d. matrix and hence |0| > 0. Also since

n > d + k, the d + k columns of A are linearly independent with probability one because

the probability of n dimensional random column vectors of y lying in any linear subspace
of R's with dimension I n 1 is zero (with respect to the Lebesgue measure on R~dit

So, A"Q-lA is a p.d. matrix and hence |A"Q-lA| > 0. Then from (4-8) it follows that

|yTQ-ly p'O2-1p| > 0.

To integrate the expression in (4-6) with respect to E, we use the following definition

of Inverse Wishart distribution.

Definition 3. [29, p. 85] Let 8 be u p x p p.d. matrix;. Then for some m > p the p x p

random matrix: W is said to have a Inverse Wishart distribution with parameters m and 8

if the p.d.f of W1 (wi!th, respectt to Lebesgue measure on, E. I, restrictedl to th~e set wrh~ere

W > 0) is given by

m+p 1-1 -1

mp p(p-1) m ]/ .-i)
2 2 x 4 102 __ -,,,

f (W, m, 8)

and we write W ~ IW,(m, 8).

Hence if a '> d + k i.e., a k > d by the above definition of Inverse Wishart

distribution, we have

2~ ) ~ x r 1/ 0( ( k 1-i)) a

whn =d + k
If a=1 d + k using (4-8i) wegl getll \

4-9)

2 4
aca+ |A hq).( -0
I i1 4d~-

JW JR

Since h(-) is a probability density function, it follows that in the case n = d + k,

which, of course, is a finite number. So, we have proved that in the particular case when

n = d + k, the posterior distribution is proper with probability one. Then, an application

of Lemma 2 of Marchev and Hobert [28] shows that for y-almost all y the posterior is

proper for n '> d + k.

We will now show that the posterior distribution is improper when n < d + k. Let

L = yTQ-ly p'O2-1p. It is easy to see that I Q-2X(XTQ-1X)-1XTQ-2 is an

idemlpotent matrix wsithl tr I 1 -X(XT-1'X)-~1XTQ- = n- k.- We also kllno that
rankly) = d with probability one. Hence, from (4-7) it follows that if a < d + k, then L

is a p.s.d. matrix i.e., there exists a vector xo(f 0) such that Exo = 0. We use this fact to

show that the integral in (4-9) diverges when n < d + k.

Since E is a symmetric matrix, the Jacobian of the transformation Z = E-l is given

by J = |Z|-(d ) [29, p. 36]. Therefore, we have

|E|exp tr(E 2)d Z exp -trZ) Z
w2 2

The matrix Z is p.d. So, by Clo d.1 -l:y decomposition, Z can be uniquely represented as

Z = LLT, where L is a lower triangular matrix (1.t.m.) with positive diagonal elements.

Th~e Ja~cobian of th~e corresponding transformationl is given by | J| = 2" nl0 1 -i+ 1[29, p.

36]. Let L = (li, 12 d). Therefore,

|f, Z| exp tr (Z) dZ= 2 | Lrn-kd-1ep -t L 1-+d
1i 1
= 2 |Lt-kd-exp -tr 1 -+L
i=1i= 1

=~~~~~ 2 x 101 1-kiL
ad~ ~ i 1 --dl~ i=1 1=

where

U = { (Un1, U21,, i2, L, Udd) E H2 ii > 0, Vi = 1, d}.

Notice that the columns of L form a basis of R d. So there exists constants bl, b2,..., bd

with b, j 0 for some i such that xo = 2, bil Suppose i'= min {ie {1, 2, .. ., d} : bi j

0}. Now consider the transformation L = (li, 12 d) O = (ox, 02, Od) Where

oi = 14 V i / i' and oi/ = zo i.e., O = LA, where

1 0 --- 0 --- 0

0 1 --- 0 --- 0

A=
0 0 --- bi/ --- 0

0 0 --- bd --- 1

Then the Jacobian of the transformation is | |A| | = | bi |-d [29, p. 36]. Note that

14, = of, -e D bo and in particular lits, = oi/i/ because oi/ = 0 for i > i'. Since

foi/ = 0 we have

d d
exp i1 01 1 -k-idL
ifi 1" i= 1

exp o s bo r- iso -d
i/i i' i>i/ i>i' i=1

exp o foT + 1c o Fo1 oi O-k-idO

vi/ ii/ i i=1

where

V = {(v1 v21, U22, 31, Udd) ER : 2 > 0, Vi / i' and A~ I > 0}.

By Fubini's theorem we can rearrange the order of integration. Notice that oi/ does not

appear in the exponential term and the only term involving oi/i/ in the above integral is

oppk-i'~. Hence the above integral dliverges since

oppk-[ i R i'~i')I(bs >o 0) +Is os,)rIcbs, 0u) does, = 0

So, we have now proved that the posterior is proper for y- almost all y if and only if
n~d+k. O

As we mentioned in the introduction that Fernandez and Steel [12] gave a proof

of propriety of the posterior density xr(P, Ely). A byproduct of our alternative proof is
a method of exact sampling from xr(P, Ely) in the particular case when n is d + k. We
describe the method now.

Let xr(q, p, Ely) be the joint posterior density of (q, p, E) given by

Then from (4-5) it is easy to see that

/( "(q, P, Ely)dq = xr( Ely). (4-11)

So iid draws from xr( Ely) can be obtained by making i.i.d draws from xr(q, P, Ery) (and

then just ignore the q component). Draws from xr(q, p, Ely) can be made sequentially via

the following decomposition

In the special case when n = d + k, from (4-10) we know that xr(qly) = nL, h(qi). So, in

this case an exact draw from xr( Ely) can be made using the following three steps:

(i) Draw 41,q2,...q,, independently where qi ~ h(qi).

(ii) Draw E ~ IWd[n k, (yTQ-ly yTQ-1X(XTQ-1X)-1XTQ-ly)- ]

(iii) Draw PT ~ N~d~k NT -1X(XTQ-1X)-1, E, (XTQ-1X)-1)

Standard statistical packages like R (R Development Core Team [36]) have functions for

generating random matrices from the Inverse Wishart distribution. One way to generate

Z ~ Nm,r(p-, E, T) is to first independently draw Ze ,pT hr steihrwo

p- for i = 1,..., m. Then take
ZT

ZT
Z= E2

ZTn

4.3 The Algorithms

In this section, we develop the DA and the Haar PX-DA algorithms for the posterior

density (4-3). We first develop the DA algorithm in Section 4.3.1 using the latent data

q = (qi, q2,., q n) and the joint posterior density of xr(q, p, Ely). We then derive the Haar

PX-DA algorithm in section 4.3.2. In the special case, when the observations, yi's, are

assumed to be from a multivariate Student's t location-scale model, van Dyk and Meng

[51] developed the marginal augmentation algorithm, which is a modified version of the

standard data augmentation algorithm, for the density (4-3). Hobert and Marchev [17]

have shown that when the group structure exploited by Liu and Wu [27] exists, marginal

augmentation algorithm (with left-Haar measure for the working prior) is exactly same as

the Liu and Wu's [1999] Haar PX-DA algorithm. In section 4.3.2 we show that a similar

group structure can be established for analyzing the posterior density in (4-3) and so

marginal augmentation algorithm is the same as the Haar PX-DA algorithm in our case.

4.3.1 Data Augmentation

We now describe the basic data augmentation (DA) algorithm. The DA algorithm

simulates a Markov chain with Markov transition density

where xr(P, E1q, y) and xr(q|4, E, y) are the conditional densities obtained from the joint

density xr(q, p, Ely).

Note that, k(P, E|4', E')xr(P', E'|y) = k(P', E'|4, E)xr(P, Ely) for all (p, E), (P', E') E

E"'' x W, i.e., k(P, E|4', E') is reversible with respect to xr(P, Ely). It then immediately

follows that xr(P, Ely) is the invariant density for the DA algorithm, i.e.,

Since the Markov transition density, k(P, E|4', E') is strictly positive (with respect to the

Lebesgue measure on E"'' x W), using similar arguments (that we used to show the AC

algorithm is Harris ergodic) as in Chapter 3, it can be shown that the DA algorithm is

Harris ergfodic. Then from the discussion in C'!s Ilter 2, it follows that the ergfodic averages

based on the DA algorithm can be used to estimate posterior expectations.

A single iteration of the data augmentation algorithm first uses the current state

(P', E') to generate q from xr(q|4', E', y) and then draws the new state (p, E) from

xr(P, Elq, y). Simulating from xr(P, clq, y) can be done sequentially by first drawing

E from xr(E|q, y) and then drawing P from x ( | E, q, y). From section 4.2, we know

that conditionally P|E, q, y follows a Matrix Normal distribution and the conditional

distribution of Elq, y is an Inverse Wishart distribution. Conditional on (P, E, y) qi's are

independent with

gi|#F, E, y in hfin x (x -y)E-(x -y

In the particular case when h(-) is Gamma (" ) ine., when yi, y2, .. "r i" i nrTOrS~umed to

be observations from Multivariate Student's t regression, we get

ind v + d v + (P Xi yi)TE-1(PTxi yi)
2' 2

4.3.2 Haar PX-DA Algorithm

In this section we derive the Haar PX-DA algorithm using the left Haar measure as

described in Hobert and Marchev's [2008] Section 4.3. From Section 4.2, we know that

"(q y) c |X Q-] X-I |y Q- y- QT- X(X Q--X- X)- _, X Q y qhi)
i= 1

As in chapter 3, let G be the multiplicative group R,+ where group composition is defined

as multiplication. Again, the left Haar measure on G is vl(dg) = dg/g where dg denotes
then Lebesgue m measure on Rm .Lt G1 act on R" through component wise multiplication,

i.e. gq = (ggi, ,.i.; ). Then it is easy to see that Lebesgue, measures on Rn is relatively

left invariant with multiplier X(g) = g". In order to construct the Haar algorithm, we need

to verify that m(q) = faxr(.i;|
given y7 can also be written as
d n-k d

(ur~uXTQ-1X yTQ-ly y Q-1X(XTQ-1X) X Q-1 y ii2 4i4 q4hB

Therefore it follows that xr(l-i.;|<) = c3 =1, */*]~ ) Where c3 1S a COnStant that does not

depend on g. Hence,

m~~~q)~ = 3O )O = C3
i= 1 i= 1

In1 order to show that n(ql) < co, sup~pose~ x ~ -h( ); x > 0 and consider th~e standard

noninformative prior, 1/o-, for the scale parameter, o-. Then the corresponding posterior

distribution is proper since

/1 x da 1 "x
-h()hyy(y=-
o o- o- o- xo

10

Then by Lemma 2 of Marchev and Hobert [28I it follows" thatI fo =1o~)$<0 for~ ~ ~ ~ ~ ~ ~ ~~~~'\ al o x,-- .> n n>1snc ,, a be viewed as the posterior density of the scale parameter, o-, when the standard noninformative prior, 1/o-, is combined with the likelihood function of a iid observations,xl, x2., *,n, from the scale family -h( ). He~nc~e, it follows that mr(q) < 00. Consider the following univariate dlensity on RW_ e,(g) = lit) m(q) i= 1 Then one iteration of the Haar PX-DA algorithm, which is a modified version of the DA algorithm, consists of the following three steps: (i) Draw q ~ xr(q|S', E', y) (ii) Draw g ~ e,(g) and set q' = gq (iii) Draw p, E ~ xr(P, E1q', y) The Markov transition density of the PX-DA algorithm can be written as where R(q, dq') is the Markov transition density induced by the step 2, that takes q q' = gq. In the special case when we have multivariate Student's t data, it is easy to see that the density e,(g) is Gamma (n f-). So, in this particular case, the step 2 of the Haar PX-DA algorithm is simply a draw from a Gamma distribution. Our results in the next section are specific to multivariate Student's t regression. 4.4 Geometric Ergodicity of the Algorithms In this section we prove the geometric ergodicity of the DA and Haar PX-DA algorithms. In C'!s Ilter 2, we mentioned that the geometric ergodicity of a Markov chain can be proved by establishing a drift condition and an associated minorization condition. In this section, we show that these conditions can be established for the DA algforithm and the Haar PX-DA algforithm. We now prove used to establish the drift condition. Lemma 1. Sup~pose that P is a p.d. matrix: and that P - vector x. Then, the following lemma that will be XXT is a p.s.d. matrix: for some X"P- x < 1. Proof. Consider the matrix Calculating the determinant of the above matrix twice using the Schur complement of P and 1, we get the following identity |P|(1 x P- X) = |P xx |, i.e., |P XXT | x P- X = 1- |P| Since P is a p.d. matrix, |P| > 0. Similarly, since P-XXT is a p.s.d. matrix, |P--xxT| > 0. Then from the above identity it follows that X"P-1X < 1. O The following lemma establishes the drift inequality for the DA algorithm. Lemma 2. Let V(P, E) = CE (yi PTXi)TE-l(yi PTX,). Then for the DA dy..~ rithm we have n+d-k n+d-k E(V(P, E)| I', E') < V(P', E') + v.. v+d-2 v+d-2 Proof. Recall that the Markov transition density of the DA algorithm is JRY ~(p, clp', C') So we have E [E {V~(P, E)|q4,y}|O'I, E', y] G(V(P, C)IP', C') (4-12) To calculate the above conditional expectations we need the corresponding conditional distributions. The required conditional distributions, as derived in the previous sections, are the following PT | q, y"1dr(-T ~, Nd E~q y~ Ie "Q-y p and ind (V+dv+(Ti yi T 1(P~ T i q4|#, E, y ~0 =12.. 2 2 Starting with the innermost expectation, we have E(V(P, E)|E,q, y) i= 1 i=1i i=: 1i i=- 1~ 1 i= 1 To/ cacuat th bv xetton eueteflo in prpryfmtixnra distribution.i Wk o -1pul T), Wh aoere by V~ W,(,9, we meahefllwng thapety V f h as x dmnsonal noncentral Wishart distribution with m degrees of freedom, covariance matrix W and with noncentrality parameter 6 [2, chap 17]. In this case, E(V) = mW + 5. If 6 = 0, we wi that V has a central Wishart distribution and we write V ~ W,(m, 9). So, we have E [E {E(V(p, E)|E, q, y)|q, y) |', E', y] . E [ E- |E, q, y] dR + pEI-1p-1 and hence, E(V(P, E)|E, q, ) =yTC yi i= 1 i= 1 Ij-~ p x4+ x[da + pE- p ] x4 i= 1 i= 1 i= 1 p 1Xi) E 1(y To calculate the second level of expectation in (4-12), we use the fact that if X~ IW,(m, 8), then X-l ~ W,(m, 8) and E(X-1) mO [29, p. 85]. Thus, E [E(V(P, E)|E,q, y)|Iq, y i= 1 p- xi)+di X 2Xi i= 1 p x4)T E(E-ly |~) (y n l xSi) + dC xf axs. i= 1 (4-13) p 1x ) (y Q- y p'O- p) (y Now, yTQ ly-p p ~T yTQ y p p S-1 2-T2pl p1 q ysy + p (X Q- X)p 2p 1X Q ly n i= 1 n i= 1 a~ ( q ys- p x )(ys i= 1 p 1X ) . (yi p 1X ) (y Qly pl p) ( p x j=1 (Yi I-1TXi) ( j= 1 4i 4i (ys p x ). T p xy(yy- pxj) Since we assume that n '> d + k, from Section 4.2, we know that yTQ-ly p'O2- p is poitve efiit matrix,,c:, withcl proablity:l c 1. So pn (y~ p xj)(yj p'xj) = yi Q-' y)-p 0- p11 1 is a p~d. matrix w-ith probability 1. Also, it is straigh~tforward to see is a positive semi definite matrix. So, an application of Lemma 1 yields (yi p Txi) (y Q-1 y p 0- )- (yi p x ) I -. (4-14) Since we assume that r(X) = k, it implies that X Q- X is a p.d. matrix and obviously xy/ x.:IT is a p~s~d matr~ix. So ano0ther7 application? of ~Lemma. 1 give~s xfaxe8 = x (X QolX- X) x j=1 < (4-15) Therefore, we get E(V( n, E)|q, y) < (n + d k) 1 i= 1 Now, recall that ind 1 V + d v + (Txi yi)TE-1(PT4 x -i ys) , q4|#, E, y ~ 0 =12.. 2 2 So, using the fact that if w: ~ Gamma(a, b) then E( ) = finally, we havei E(V(P, E)|#', E') 5(s-# i E)-(s-# i i= 1 which proves the lemma. O The following lemma establishes an associated minorization condition. v+ (p' x, yi) (E')- ( Xi -- yi)\ v + (P' xi yi) (E')- ( X, yi) Lemma 3. Fix: 1 > 0 and let S {((, E) : V(P, E) < 1}. Then the M~arkov transition 1, .:1;i of the DA rlly.>rithm k(P, E|4', E') .el.:. the following minorization condition k(P, E|4', E') > ed(P, E) V(P', E') ES where the 7. ;. -.1 / d (P, E) is given by at~ So gtd9(i ) d r and e = ( fo g(t)dt)". The function g(-) is given by i/t (v+d v vd v+1 9() F 2 2' (,* whAere q* = iSlogl 1+ andr I(a~blx) denotes the Gamlrmajab) I r.' .1/ evaluated at the point x. Proof. For i = 1, 2,. ., n, define p'x,) (E)- (ys PTXe) < 1}~ Si = {((, E) : (ys Clearly, S C Si for i 1, 2,. ., n. Recall that, xr(q|0, E, y) is product of a Gamma densities. So, Yi ) (E')- ( X, mnf r(q|S',E', y) (P',C')ES Then by Hobert's [2001] Lemma 1, it straightforwardly follows that (q|O', E', y) > gii)ca Vci', E') t S. i= 1 V + (P' Xi .P,')~= v + d i= 1 i=1 (P',C')ES 2 ' Hence, the proof follows because Fr-om results stated in chapter 2 it then follows that the DA algorithm is geometrically ergodic as long as Ithe coefficient of V(B', E') in Lemmnra 2, i.e., "- is strictlly less than 1. We state this in the following theorem. Theorem 4. The DA ,Ily.>rithmn is ~i..- tre. I :. .rli ergodic if 0 < < 1 i.ec., a2 < v+k-2. Hobert and Marchev's [2008] Proposition 6 shows that the Haar PX-DA algorithm is at least as efficient (in efficiency ordering) as the DA algorithm. Using similar arguments as in Corollary 1, we can show that geometric ergodicity of the DA algorithm implies that of the Haar PX-DA algorithm. Hence we have the following corollary. Corollary 2. The Haar PX-DA rlly. rithm is at least as efficient as the DA rlly.>rithm. Also, the Haar PX-DA rlly.>rithm is geomen,... ril;i ergodic if a < v + k 2. Remark 4. Our result in Lemma 2 matches with M~archev and Hobert's [2004]l Theorem 1' in the case when k = 1. We also can prove the geometric ergodicity of the Haar PX-DA algorithm by directly establishing a drift and minorization condition for it. We actually can use the same drift function V(P, E) = CE (yi PTXi)TE-l(yi PTX,) to establish a drift condition for the Haar PX-DA algorithm. Lemma 4. Let V(P, E) = CE =(yi PTXi)TE-l(yi PTX,). Then for the Haar PX-DA rlhi ~.:athm we have (n +d k)(v +d)(n 1) av(n +d k) n(n 1)(n + d k)v(v + d) E ( V (, E) | ', E')I (nv 2)(v +d 2) nv 2 (nv 2)(v + d 2) Proof. Recall that the Markov transition density of the Haar PX-DA algorithm is - yi) E-1(P X, - So we have V(#, E)k(p, E|4', E')dpdE JW JR E [E { V(#, E)|q, y) | ', E', y] E(V(P, E)|# ', E') From Section 4.3.2 we know that q' gq where giq ~ Gamma (n ~-). We substitute q by q' in (4-13) and then straightforward algebra shows that (vi p xi) + xOx i= 1 n k) 9 p 1Xi) (y Qly pl p) (y E(V(P, E)|q1, y) E(V(P, E)|q, y) = (n k) (ys l i = 1 i i= 1 p'1Xi) (YTQ- y p'1Xi) (YTQ- y p'O- p)- (y p'O- p)- (y p- xs) + di x Ox i= 1 From (4-14) and (4-15) it then follows that v~n d -k) "1 nu 2 q4 i= 1 E(V(P, E)|q, y) v(n +d k) g av 2 q::4j~ i= 1 j i Since conditional on (P, E, y), qi's are independently distributed with in s v +d V + (PTxi q4|#, E, y ~ 1,2,...,n, we have E(V(P, E)|4', E') v(nn +- d k) n+ gy j~~--> i= 1 j i v(n + d k) [ v + -vd V + (P'TXi - av-2 + 2 i=1si v~ V+ (Pl'z Yi )T(E')-1(P'TXs -yj)TE'- (P'Txy -yj>Y) E q, y yi>)2 v(n + d k) I (n 1) (v + d) (' )CI-PTX nu 2 v(V+ d 2) - which proves the lemma. O Our minorization condition for the DA algorithm straightforwardly generalizes to a minorization condition for the PX-DA algorithm. Lemma 5. The M~arkov transition I7.* .:1/ of the Haar PX-DA rlly..,7:thm k(P, E|4', ') al. J. the following minorization condition where the 7. ,: -.1 / d (P, E) is given by and e, S and g () are as .I~ .,: .1 in lemma S. Together, Lemmas 4 and 5 prove the following theorem. Theorem 5.'lt The Har PX-DA rly. rithmr is geomenl~:...rl c/ ergodi if 0 < (n~'d-k vd n1 As a corollary of Theorem 4 (see Corollary 2) we know that the PX-DA algorithm is geomnetric~ally erg~odic if < 1. At first it might appear that Theoremr 5 is a better result than Corollary 2. But, we now show that it can never happen that both of the followingf inequalities hold together n + d-k (n d -k)(v d)(n -1) > 1 < 1. (4-16) v+ d- (nv -2)(v d -2) N\ote that (n d-)(v'd)(n-) n '+d-k if and only if (" ")( < 1 i.e. nr < 1 + "-2. So, if (nv- 2)(v+d- 2) v+d-2 (nv-2) d (4-16) holds then the following must holds n+d-k v-2 > 1; ,a< 1 (4-17) v+d-2 d which, in particular, we~ that v has to be '?N r;i than 2. Now (4-17) holds if and only if v-2 n~v+k-2andn<1+ But, clearly, the above two inequalities can't hold together. We, of course, would like to be able to wi that the DA algorithm ( or the Haar PX-DA algorithm) is geometrically ergodic for all (d, k, v, n)-tuples such that d > 1, k > 1, v > 0 and a > d + k. Note that the conditions in Theorem 4 and Theorem 5 are upshots of the particular drift function V(P, E) and the inequalities (4-14) and (4-15) that we used to prove the drift conditions. In our opinion, to make substantial improvement of Theorem 4 and Theorem 5, we either have to consider a different drift function or resort to a altogether different technique of proving geometric ergodicity other than establishing drift and minorization condition. Either of these two would require us to start from scratch. CHAPTER 5 SPECTRAL THEOREM AND ORDERING OF MARK(OV CHAINS This chapter is divided into two sections. In the first section we present a brief account of operator theory on Hilbert space. In particular, we define the spectral theorem for a bounded normal linear operator on Hilbert space. In the second section we show how these results from functional analysis are used in studying the theory of Markov chains. 5.1 Spectral Theory for Normal Operators We begin with the definition of a Hilbert space. Definition 4. A complex vector space H is called pair of vectors y and b in H is associatedd a complex number (g, h). called the inner product of g and b. .such that the followings hold: (g, h) =(h, g) (g + h, k) = (g, k) + (h, k) if h, k E H (a~g, h) = co(g, h) if 9, h EH and n~EC@ (h, h) > 0 for <<11 h EH and (h, h) = 0 only if h = 0. It can be .;-.le;, shown that the above inner product induces a norm on H I. I;,.. l by | |h | = (h, h) E If the r. tel. e:11 normed space is complete. it is called a Hilbert space. Let H he a Hilbert space over C and T : H H he a linear transformation, called a linear operator. The operator T is said to be bounded if there exists an At > 0 such that ||Th|| < Af||h|| for all h e H, where the norm, || ||, is as defined above. Let B(H) be the collection of all linear, bounded operators from H into H. For Te E (H), define the operator norm of T by i ~i||Th||ll htlhO It is easy to see that ||T|| = sup{||Th|| : h e H, ||h|| < 1} = sup{||Th|| : h e H, ||h|| = 1}. It can he shown that B(H) with the above norm is a Banach space (i.e., it is complete). If S, Te E (H), the composite operator ST E B(H) is defined by (ST)(h) = S(T(h)) h e H. In particular, powers of Te E (H) can he defined as To = I, the identity operator on H and T's = T,-1, for n = 1, 2, 3, .. We can easily verify the inequality || ST|| < || S|| || T|| Now, we define the adjoint of an operator. Let Te E (H). Then using Riesz representation theorem it can he shown that there exists a unique T* E B(H) for which (Ty, h) = (g, T*h) for all g, h e H. The operator T* is called the adjoint of T [47, p. 311]. Some properties of adjoint operators are listed below: ||T*|| = ||T|| T** =T ||T*T|| = ||T||2 and if S also belongs to B(H) then (ST)* = T*S*. We now can define different types of operators. An operator Te E (H) is said to be normal if TT* = T*T, self-adjoint if T* = T, unitary if T*T = I = TT*, projection if T2 =T Clearly, self-adjoint and unitary operators are normal. Note that, if T is self-adjoint then (Th, h) = (h, Th) for all h e H. But, from Definition 4 we know that (h, Th) = (Th, h). So, for self-adjoint operator T, (Th, h) is real for all h e H. We now state a useful uniqueness theorem [47, p.310]. Theorem 6. If Te E (H) and if (Th, h) = 0 for every h E H. then T is the zero-operator. Remark 5. The above theorem would fail if the .scatler field were RW instead of C. For exravple. consider the mardrixr T = < t s As a corollary of the above theorem we get an alternative definition for a normal operator. Corollary 3. An operator T on H is a normal operator if and only if ||T*h|| = ||Th|| for every h E H. Proof. Suppose that || T*h || = || Th || for every h E H. Then || T*h1 || 2 h12 and (T*h, T*h) = (Th, Th). This yields (TT*h, h) = (T*Th, h) i.e., ((TT* T*T)h, h) all b in H. By the above theorem we get TT* = T*T, i.e., T is normal. Conversely, if T is known to be normal then we can write ((TT* T*T)h, h) all b and now following the above reasoning in the reverse direction we get ||T*h|| for every h E H. so = for 0 for = |Th|| O Suppose Te E (H). The operator T is invertible if there is an S E B(H) such that TS = I = ST. The operator S is called the inverse of T and we write S = T l. Suppose R(T) denotes the range of T, i.e., R(T) = { Th : h e H}. Then T is invertible if and only if R(T) = H and T is one-to-one. The spectrum, o-(T), of the operator Te E (H) is defined as follows o-(T) = {A : T AI is not invertible}. Thus Ae o (T) if and only if at least one of the following two statements is true (i) The range of T AI is not all of H i.e., T AI is not onto. (ii) The operator T AI is not one-to-one. If (ii) holds, A is said to be an 1.:I ,: al;,.- of T and in this case there exists & / 0 such that Th = Ah. We call h an eigenvector corresponding to the eigenvalue A. It is difficult to make a general statement about the spectrum of an operator. We define compact operator later in this section. We will see that we can ;?, quite a lot about the spectrum of a compact operator. The complement of o-(T) is called the resolvent set of T and is denoted by p(T), i.e., p(T) = C \o-(T). The proof of the following theorem can be found in text books on functional analysis [see e.g. 9, p. 83]. Theorem 7. If ||T || < 1, then (I T) is an invertible operator. Furthermore, n= o Corollary 4. Let A be a complex; number such that ||T || < |A|. Then from the above theorem it follows that T AI is an invertible operator and its inverse is noi Corollary 5. Suppose that T is an invertible operator and S is another bounded linear operator such that ||T S|| < ||T-1|| ' then S is an invertible operator. Proof. We have || T- S I = IT-(T S)I < IT- IIT SI < 1 by our hypothesis. Thus, by above theorem, I (I T-1S) = T-1S is an invertible operator. Let us denote the operator T-1S by D. Then S = TD and so D- T-l is the inverse of S showing that S is an invertible operator. O Remark 6. Cor .ll.r,, a tells us that a(T) is a bounded set for r,.;; bounded linear operator T because a(T) C {A E C: |A| < ||T||}. Theorem 8. For r ;, bounded linear operator T the set a(T) is a closed, bounded subset of C. Proof. We already know that a(T) is a bounded set. We only need to show that a(T) is a closed set. We will show that p(T) is an open set. Let A be any point in p(T), i.e, (T AI)-l is a bounded linear operator. Choose p such that |p 1- A| < ||(T AI)- ||-l and let T' = T AI, T" = T pl. Now T' is invertible and From Corollary 2, it follows that T" is an invertible operator. Thus, p E p(T) and in particular it follows that i.e., any A s p(T) has a neighborhood contained in p(T). Hence, p(T) is an open set. O It can also be proved that for Te E (H), H / {0}, the spectrum of T, o-(T), is not empty [47, p.253]. The spectral radius, rT, of T is defined as rr, = sup{|A| : A E o-(T)}. So, rT is the radius of the smallest disc (with center at 0) containing the spectrum of T. From Remark 6, we see that rT < ||T||. We now state the spectral radius theorem. Theorem 9. (Sp~ectral Radius Theorem) If Te E(H). the spectral radius. rT. of T r, = lim || T'sl || in f ||T'j The operator norm of a normal operator is same as its spectral radius. We state this as the following theorem. Theorem 10. Let Te E (H) be a normal operator. Then rT = ||T||. If T is self-adjoint then a(T) is a non-empty compact subset of RW and is therefore contained in some smallest closed interval, [nz(T), Af(T)], of RW. Also, ni(T) and Af(T) can he expressed in terms of the inner product on the Hilbert space. These along with some other properties of self-adjoint operators are stated in the following theorem Theorem 11. Let Te E (H) be closed interwel containing a(T). Then (i) nz(T) = inf (Th, h) (ii) Af (T) = sup (Th, h) and (iii) ||T|| = .sup |(Th, h)| = mesxr(|m(T)|, |AF(T)|) We define an operator, T, to be positive if (Th, h) > 0 for all h eH and we write T > 0. An operator T is positive if and only if T is self adjoint and a(T) c [0, 00) [47, p. :330]. It is easy to see that both TT* and T*T are positive operators for any Te E (H). If S, Te E (H) are two self-adjoint operators, we ;?-, S > T if and only if S T > 0. If we consider the scalar field to be RW instead of C, there may exist a positive operator T that is not self-adjoint. The following theorem can he found in [47] p. :331 [also see :37]. Theorem 12. Every positive operator T E B(H) has a unique positive square root S .such that S2 = T (184 U)C Urite S = Ti. If T is invertible. .so i~s S. Now, we describe the spectral theorem for normal operators which defines the operator f(T) where f is a bounded Borel measurable function defined on the spectrum, o-(T), of a bounded normal operator T. A nice exposition on spectral theory is given in Devito [9]. He starts with deriving the spectral theorem for operators on finite dimensional vector space in chapter 5 and ends in chapter 7 with spectral theorem for unbounded self-adjoint operator. There are various approaches to derive the spectral theorem for a bounded, normal operator. In this chapter we define spectral theorem following Rudin [47]. First, we define resolution of the .:ll. ,1.:1; Definition 5. Let B(R) be 8 : B(n) i (H) with the following properties: (i) (4) = 0. the zero operator. 8 (R) = I (ii) each 8(w) is (iii) S( (n w') = S (w) S(w') (iv) If w n w' = 4 then S( (U e') = S (w) + S (w') (v) for every g, h E H. the .set function. Seas. I' It...lI by is a complex measure on B(R). Since each S(w) is a self-adjoint projection, we have Simz( o) = (S(Lo)h, h) = ||(W~h||2 for all h E H. So from (v) it follows that Simit is a positive measure on B(R) for each hE H. Theorem 13. (Sp~ectral Theorem) Let T ES (H) be a normal operator. Then there exists a unique resolution of the .:1. ,:./:7;, 8, on B(a(T)) such that (Tg, h) =m AS,,(dA) for all g,h h H L/et f be a bounded Borel function on a(T). The operator f(T) is I. I;,.. l to be the unique operator ;oril r;.,ii l The above theorem justifies the following notations and f (T) = f I(A)Sl~dA) We refer to the above 8 as the spectral decomposition of T and we denote it by Sr. We now give an example. We know that any linear operator, T, on a finite dimensional vector space, V, can be represented as a matrix [9, p. 157]. We define T to be a diagonalizable operator if its matrix representation is a diagonalizable matrix. From [9, p. 167] we know that if T is a normal operator on V, then T is diagfonalizable. In particular, if we take V to be RW", then T can be represented by n x n diagonalizable matrices. Assume that all the eigenvalues of T are distinct and let Ai, A2, *, abe the distinct eigenvalues of T with corresponding eigenvectors P1, P2,..., P,. Hence, a(T) = {Az, a2, *, n} and B(a(T)) is the powe~rset, of e(T(.). For Ae B t(c(T)), we define Sr(Al) = i:AgEA Piri'. If A = weV define Sr(#) = 0, the zero operator. Since the eigenvectors Pi's are mutually orthogonal [9, p. 157], it is easy to see that Er(-) satisfies all the properties stated in Definition 5. Also note that, i= 1 1\ore generally, if T is a normal operator on a finite dimensional Hilbert space then Sr(A) is the orthogonal projection onto union of all eigenspaces for {Ag}, where As'sS are eigenvalues of T that are contained in ,4. The operator f(T) in Theorem 13 is of course a bounded operator and it can he proved that || f(T)|| < sup{| f ()| : Ae E (T)}. We also have the following theorem which determines the spectrum of the operator f(T). Theorem 14. Let Te E (H) be a normal operator and let f : a(T) C be continuous. Then f (T)* = f(T) || f(T)|| = mp{| f (A)| : Ae E (T)}. As a corollary of the above theorem we get the following useful result. Corollary 6. Let Te E (H) be a normal operator. Then |I ||= ||T ||' for every positive integer n. Now we introduce an important class of linear operators called compact operators. Recall that a subset ,4 of a Banach space is compact if every sequence {a,z} in 4 has a subsequence converging to some element of A. Definition 6. We n;, Te E (H) is compact ifT maps the unit ball. {h eH : ||h|| < 1} into a .set in H whose closure is compact. Retherford [37] gives a nice, concise description of compact operators. For a compact operator T, the spectrum, a(T), is at most countable, and has at most one limit point, namely, 0. Also, any non-zero value in the spectrum is necessarily an eigenvalue of T. We state these results in the following theorem ([8] p. 214, [9] chapter 5 ). Theorem 15. Let Te E (H) be a non-zero, compact, self adjoint operator. Let {A,} be the distinct, non-zero .:I c. ,:;rle;, of T arr rtg J:I~~ in such a way that I11 > I12 > i jfOri )j. and let {h,} be the corresponding sequence of orthonormal eigenvectors. If the sequence does not terminate, then lim |A,| = 0. Moreover, for h E H, ntOO n= 1 L/et G, be the eigenspace of As, i.e., G, is the closed linear span of all the As 's in {hm} having As as ~. .9. :;...le;, Then each G, has finite dimension and its dimension is the number of times An is repeated. F.:,..elle; the spectrum of T, o-(T), is ,.;.-l;; {A,} together with {0}. Some useful properties of compact operators are listed below: T is compact if and only if T* is compact. If T is compact and S E B(H) then TS and ST are compact operators. T is compact if and only if T*T is compact. T is compact if and only if TT* is compact. Every linear operator on a finite dimensional vector space is a compact operator and its spectrum coincides with the set of eigfenvalues of the operator. Another example of a compact operator is the Hilbert-Schmidt integral operator ([47] p. 112, [8] p. 267 ). Let (R, B(R), p) be a measure space. Suppose, t : O x R C be such that |tI (x, y) |2iilly< L Then the associated Hilbert-Schmidt integral operator is the operator T : L2( @ L2(2; 0,) giVen by (T f )(x) = tx )f()d~) 5.2 Application of Spectral Theory to Markov Chains Recall from C'!s Ilter 2 that P(:r, dy) is a Markov transition function on (X, B(X)) with invariant probability, measure, x Lt L (x) be the vector space of real-valued, measurable, mean-zero functions on X that are square-integrable with respect to xr, i.e., Ifw defin an;, inne produc~,,t on~2, L, (x y (,g) = x f(:r)g(:r)iT(dX) then L2, (x i Hilbert space with its norm given by || f|| = ( f, f) The Markov transition function, P(:r, dy) defines an operator, P, on L (xr). The operator P takes each f eL (xr) to the function, P~l fa\ (eL (x)), deindasfllows, By Jensen's inequality it follows that ||P|| < 1. So P is a bounded linear operator on L (x). The vector space L2- (x) is a subspa Ce of2(x) and' L (x) is the- spac tha~t is orthogonal to constant functions. We later describe the reason why we consider P to be an operator only on L2, (x)intead of the whole,, spac 2(). We would also like to mention that, unlike the previous section, here we consider the scalar field to be RW because in statistics we are mostly interested in real valued functions. In this case, the inner product is symmetric in its arguments i.e., ( f, g) = (g, f). Fr-om C'!s Ilter 2, we know that the Markov chain A (and so the transition function P(:r, dy)) is reversible if and only if for all f, y E L2(;~ (P f, g) =( f, Pg). It is easy to show that the above definition of reversibility is unchanged if the space L2(,) is replaced with L (xr). So, the Markov chain A is reversible if and only if P is a self-adjoint operator on L (xr). The operator I P, known as the Laplacian operator corresponding to P, pl .i-~ an important role in analyzing Markov chains and we denote this operator by lp. Note that the operator 1p is not one-to-one if there exists f (f 0) E L (x) such that Ip f= 0. But, from ?1. i-n and Tweedie's [1993] Proposition 17.4.1 it follows that for Harris ergodic Markov chains the only functions in L2 K) SatiSfying the equation Ip f = 0 are xr-almost everywhere constant. Since L (xr) does not include functions which are xr-almost everywhere constant (non-zero) and since we consider Harris ergodic Markov chains, the operator Ip is one-to-one in our case. But, Ip might not be an onto operator and so it might not be invertible. From the definition of spectrum of an operator we know that lp is invertible iff 1 / a(P). Since a(P) is a closed set it follows that 1 tf a(P) if and only if supa(P) < 1. Hence the Laplacian operator, lp, is invertible iff supa(P) < 1. Recall from C'!s Ilter 2 that we ;?- a CLT holds for fm (or simply f) if there exists a ,2 E (0, 00) such that, as m 00o, ~(~f Ex f) iN~(0, a"2 It is proved in [13] that a CLT holds for any function f lying in R(lp), the range of lp, and the corresponding .I-i-uspind ic~ variance in the CLT is given by v(f, P) = ||g||2_- l 12, where f = 1pg g EL (jT) So, if Ip is invertible, then CLT holds for every function feL1, (x. ls, from Teore ?,, we know that if ||P|| < 1, then Ip is alrv-ove invertible. Hence, we have the following proposition. Proposition 2. Let P be the M~arkov transition function of a Harris ergodic M~arkov chain with invariant I1:. e. -1 / If ||P || < 1, then CLT holds for,,,, evey untion fe (x .2, Remark 7. Note that, for a Harris ergodic M~arkov chain on a finite state space, lp is rl;, ate, invertible. In this case a(P) contains only the .:II no.;rle;,. of P and 1 tf a(P) because 1 is not an .:II ,.:.;rle;,. of P since we consider P as an operator on Li(xr). So, for a Harris ergodic M~arkov chain on a finite state space, CLT holds for every function fe U\L (x). Since Ip is one-to-one, even in the case when Ip is not onto we still can define the inverse operator ly] if the domain of ly D(lp ), ;?--, is restricted to R(lp) [32]. Fr-om [13] it then follows that CLT holds for all f e D(lfl) and by lemma 3.2 of Mira and G. Or-l [32] we have v(f P) = (f [21p' I]f ), V fe D (lf ). In C'!s Ilter 2, we defined the efficiency ordering due to [32]. The problem with efficiency ordering is that even when we conclude P FE Q, it might happen that v(f, P) = v(f, Q) for every f E L2(;T) that are of interest. In this chapter, we discuss conditions that guarantee v( f D,,. ~ P) < c v( f Q)f r evey eL(x (5-1) for some constant 0 < c < 1 and in this case we ;?-- P is more efficient than Q with uniform relative efficiency c. Also, we ;?i- that P is strictly more efficient than Q if vfC \, P) < v,,f,, Q)fr vry feL () (5-2) Note that (5-1) implies (5-2). But, (5-2) does not necessarily imply (5-1). Let lp and 1Q be the Laplacian operators corresponding to the Markov operators P and Q. If we assume that both lp and 1o are invertible then it follows that D(lp )= L~ (xr) D(lf ). Hence, we have v ( f, P) = ( f, [21-1 I] f), V fe ( x)- and v (j f ) = ( f [21,l I]) fo ), V fL (xr). Proposition 3. Assume that both lp and lq are invertible. Then v( f, P) < c v( f, Q) for every f E L2 K7) ;T C -1_ --1 + I ~~> 0 Proof. The proof follows from following equivalences v(f, P) m (f ,[2) -I]f ) < c(f, [21l If ) V fe L2~ + (f( c21,l I]- 21,1 I) f )>O Vfe L2r + c [21, I]-21,1 I> 0 + 2c ly1 21, + (1 c)I > 0 -1(1 c) W c 9 p + I '> 0. The above proposition is similar to Corollary 5.1 of Mira and G. Or-l [32]. Remark 8. In Proposition 8, if we replace the inequalities by strict inequalities and c by 1, we get conditions for strict e~ff. 7.: ,: ; For the rest of this section we assume that we have reversible Markov chains. We know that the Markov transition operator, P, of a reversible Markov chain is self-adjoint. So, we can use spectral theory to analyze reversible Markov chains. Let v( f, P) =,, 1 AS,7f~f(dA), (5-3) where 8(-) is the spectral decomposition of P. K~ipnis and Varadhan [24] show that if v( f, P) in (5-3) is finite, then there is a CLT for f with the .mon i nd il~lc variance, v( f, P) given by (5-3). Now onwards, we denote 87,7(-) by Syp(-). Recall from Chapter 2 that P is said to be better than Q in the efficiency ordering written P FE Q, if U f, P) I U f, Q) for every f E L2(x). It can be easily shown that the definition of efficiency ordering is unchanged if the space L2 Kr) ;,,1S, repace wth L (x). It is tempting to conclude, from (53) and th~e spectral thleoremn that P >E Q Iiff > ~. But, siinc~e o-(P) c [-1, 1], the function h(x) = is not necessarily bounded on o-(P) (o-(P) C [-1, 1] because ||P|| < 1). Hence we might not be able to use the spectral theorem to define h(P). The following definition is from Mira and G. Or-l [32]. Definition 7. Suppose, P and Q are two M~arkov transition functions with invariant I4~ll.:l..1l...;i measure xr. Then P is said to be better than Q2 in the covariance ordering written P >> Q, if Q > P on L (r) . It is easy to see that any Markov operator, P, dominates I, the identity operator, in Peskun ordering (Peskun [35], Tierney [50]). Also, Tierney [50] shows that Peskun ordering implies covariance ordering. So, P >> I, i.e., lp = I P is a positive operator. From our discussion in the previous section it then follows that there exists a unique square root 1) of Ip. K~ipnis and Varadhan [24] show that v(f, P) < 00 if and only if f E R(lf). Note that, this does not contradict the result of Gordin and Lif~sic [13] since D/l_) C R~lf). Mira, and/,, Geyrl [32]prov that P Fi Q if and only if P FE Q. We give an alternative proof for P >> Q implies P FE Q. The following theorem is from Bendat and Sherman [3] [see also 32, Theorem 4.1]. Theorem 16. Let h(x) be a bounded Borel measurable function on (lI, I2). A necessary and sufficient condition for h to have the 1/* */'' l ;, that h(S) < h(T) for all bounded, self-adjoint operators S, T with S < T and o-(S),o-(T) c (II, I2) iS tati h iS t'il: O (II, I2) Gnd Can it,'.ibil.:. .ill/ be continued into the whole upper-half plane, and represents there an r,:...;let..- function whose ie.:lrl.:..rry part is non-negative. Bendat and Sherman [3] mention that the function ax + b h(x) = with ad bc > 0 cx +d satisfies the condition of Theorem 16 either in x > -d or x < -. We now prove the following theorem. Theorem 17. Suppose that P and Q are two reversible M~arkov transition functions. If P Fi Qc then P FE -c) Proof. For E (0, 1) define 1 + x Comparing with the function h(:r) = @ ,w ae db 2( )an >1 o for e E (0, 1) the function h,(:r) is analytic on a(P) and a(Q) and it satisfies the conditions of Theorem 16. Since, we are assuming that P 1 IQ i.e., Q '> P, by Theorem 16 we have h,(Q) > h,(P) for ee (0, 1). For Ee (0, 1), h,(:r) is also a bounded Borel function on [-1,1]. So we can use the spectral theorem to define the operator h,(P). By spectral theorem, we have (h,(P) f, f) =~~ he(A)Syp~A). Let a+(P) = a(P) 0 [0, 1] and a-(P) = a(P) 0 [-1, 0) then we can write Note that, d 2A For Ae e*a(P), h,(A) > 0 and he(A) ? as e 1 Hence, by 1\lonotone convergence theorem, we have h.(A)Syp(dh) T1 1 /l + (P) + (P) For Ae t (P), he,(A) i \as e 1 and 0 < he,(A) < 1. Since S, 18yp(dAX) < || f ||2< 00 by Dominated convergence theorem we get, /,~ /,"EF~l~~ 1+Al- So, SJ,j h,(A)Syp(dAX) .Jb Sy~A asFI~, e 1. For e e (0, 1), we know that h,(Q) > h,(P) i.e. (ht(Q) f, f) >(h,(P) f, f) Vf fL (xr) 1.e., Then taking the limit e 1 on both sides we get i, 1+A_ P 1+A i \e., v Gf, Q) > vfP) Vf E O\L (x. O Proposition 4. Assume, sup?(a(P)) < 1 and sup?(a(Q)) < 1. Then P >> Q iff P FE -. Proof. It is assumed that sup(a(P)) < 1 and sup(a(Q)) < 1. So the function h.(x) = M is bounded on a(P) and a(Q). Since a(P), e(Q) C (-2, 1) and h(x) is analytic in x < 1, by Theorem 16 we have I +Qd I +P 2> Po > I I-Qd I- P Since a(P) < 1, the function h(x) is continuous on a(P). By Theorem 14, we know that h(P) is a bounded operator. Since we consider the scalar field to be RW, by Theorem 14, wve also know that h(P) is self-adljoint. Let g(x) = Then g(x) is analytic when x > -1. Let max(sup a(P), sup a(Q)) = e'. Then both a(h(P)) and a(h(Q)) are subsets of [0, ( )]. So, we can apply Theorem 16 on g to show I +Qd I +P > 4Q~2> P. I -Q d I- P Hence, we get I+Q I+P S> Po > ~ r -Q r -P The function h(x) = M is bounded on a(P) and a(Q). So, we can use spectral theorem to define the operator h.(P) = Then, byv (5 3) we have I+Q I+P > -- aP FE - I-Q r -P Hence , I +Qd I +P P ~ ~ ~ I-G)~ c ~>~P G) I d I - Now we discuss how theory of compact operators can he used to analyze Markov chains. Let P he the Markov transition function of a Harris ergodic Markov chain. Fr-on Proposition 2 we know that if ||P|| < 1, then CLT holds for every square integrable function f. It is easy to see that if P is a compact operator then ||P|| < 1. So, in order to establish CLT for a Markov chain we can show that the corresponding Markov transition operator, P, is compact. In Section 5.1, we mentioned that any Hilbert Schmidt operator is compact. Recall front ChI Ilter 2 that if a reversible Markov chain is geometrically ergodic, then the CLT holds for every square integrable function f [41]. Also we mentioned before that a reversible Markov chain, P, is geometrically ergodic if and only if ||P|| < 1 [41, 44]. Schervish and Carlin [48] and Liu et al. [26]) have proved geometric ergodicity of certain Markov chains by establishing that the corresponding Markov operator, P, is Hilbert Schmidt. APPENDIX: CHEN AND SHAO'S CONDITIONS Here we state C'I. n! and Shao's [2000] necessary and sufficient conditions for city) < 00 as well as a simple method for checking these conditions. Let X denote the nx p matrix whose ith row is XT and let W denote an n x p matrix whose ith row is WT, where Xi if yi = 0 S-Xi if yi = 1. Proposition 5. [6] The function city) is finite if and only if (i) the design matrix: X has full column rank, and (ii) there exists a vector a = (al,..., a,) with strictly positive components such that W~a = 0 . Assuming that X has full column rank, the second condition of Proposition 5 can be straightforwardly checked with a simple linear program implementable in the R programming language [18] using the -p!'!.I::") function from the "boot" library. Let 1 and J denote a column vector and a matrix of 1s, respectively. The linear program calls for maximizing 1Ta subject to W~a = 0 (J I)a < 1 (element-wise) ai > 0 for i = 1,...,n This is ah-liws feasible (e.g., take a to be a vector of zeros). If the maximizer, call it a*, is such that a,* > 0, for all i = 1,..., n, then the second condition of Proposition 5 is satisfied and cz (y) < 00. Moreover, it is straightforward to show that if a* contains one or more zeros, then there does not exist an a with all positive elements such that W~a = 0, so c (y) = 00. REFERENCES [1] Albert, J. H. and Chib, S. (199:3), "B li-o Io analysis of binary and polychotomous response data," Journal of the American Statistical A~ssociation, 88, 669-679. [2] Arnold, S. F. (1981), The Theory of Litear M~odels and M~ultimerieste A,...el;;-.: Wiley, New York. [:3] Bendat, J. and Sherman, S. (1955), "Monotone and convex operator functions," Transactions of the American Afrathematical S8... .. I ;, 79, 58-71. [4] Casella, G. and George, E. (1992), "Explaining the Gihhs sampler," The American Statistician, 46, 167-174. [5] Chan, K(. S. and G. o c., C. J. (1994), "Discussion of il ..Inl~v chains for Exploring Posterior Distributions"," The Annedls of Shetistic~s, 22, 1747-1757. [6] Chen, M.-H. and Shao, Q.-M. (2000), "Propriety of posterior distribution for dichotomous quantal response models," Proceedings of the American Afrathemati- crel 8... .:. I;, 129, 29:3-302. [7] Chib, S. and Greenberg, E. (1995), "Understanding the Metropolis-Hastings algorithm," The American Shetistician, 49, :327-:335. [8] Conway, J. B. (1990), A C'ourse in Functional Air;,;-.- Springer-Verlag, New York, 2nd ed. [9] Devito, C. L. (1990), Functional A,:tale;,.: and Litear Op~erator The .-; Addison-Wesley Publishing Company. [10] Eaton, 31. L. (1989), Group? Invaricence Applications in Stratistic~s, Hayward, California and Alexandria, Virginia: Institute of Mathematical Statistics and the American Statistical Association. [11] Feller, W. (1968), An Introduction to P,~rol~.:l..7.;i Theory and its Applications, vol. I, New York: John Wiley & Sons, :$rd ed.

[12] Fernandez, C. and Steel, 31. F. J. (1999), \!ll1 variate Student-t regression models:
Pitfalls and inference," Biometrikes, 86, 15:3167.

[1:3] Gordin, 31. I. and Lifsic, B. A. (1978), "The central limit theorem for stationary
Markov processes," Soviet Iabthematic~s. D -11.; Jr;, 19, :392-394.

[14] Hohert, J. P. (2001), "Discussion of "The art of data augmentation" by D.A. van Dyk
and X.-L. Meng," Joureml of C'ompubstional and Grap~hical Stratistic~s, 10, 59-68.

[15] Hohert, J. P. and Geyer, C. J. (1998), "Geometric Ergodicity of Gihhs and Block
Gihhs Samplers for a Hierarchical Random Effects Model," Journal of Afultivaricate
A,...el;;-.: 67, 414-430.

[16] Hobert, J. P., Jones, G. L., Presnell, B., and Rosenthal, J. S. (2002), "On the
applicability of regenerative simulation in Markov chain Monte Carlo," Biometrika,
89, 731-743.

[17] Hobert, J. P. and Marchev, D. (2008), "A theoretical comparison of the data
augmentation, marginal augmentation and PX-DA algorithms," The Annals of
Statistics, 36, 532-554.

[18] Ihaka, R. and Gentleman, R. (1996), "R: A Language for Data Analysis and
Graphics," Journal of Comp~utational and Grap~hical Statistics, 5, 299-314.

[19] Johnson, N. L. and K~otz, S. (1970), Continuous Univariate Distributions-1, John
Wiley & Sons.

[20] Jones, G. L. (2004), "On the Markov chain central limit theorem," Pir..ort.7:/i;,
Surveys, 1, 299-320.

[21] Jones, G. L., Haran, M., Caffo, B. S., and Neath, R. (2006), "Fixed-width output
analysis for Markov chain Monte Carlo," Journal of the American Statistical Associa-
tion, 101, 1537-1547.

[22] Jones, G. L. and Hobert, J. P. (2001), "Honest exploration of intractable probability
distributions via Markov chain Monte Carlo," Statistical Science, 16, 312-34.

[23] -(2004), "Sufficient burn-in for Gibbs samplers for a hierarchical random effects
model," The Annals of Statistics, 32, 784-817.

[24] K~ipnis, C. and Varadhan, S. R. S. (1986), "Central limit theorem for additive
functionals of reversible Markov processes and applications to simple exclusions,"
Communications in M~athematical Ph ;; .~ ; 104, 1-19.

[25] Liu, J. S., Wong, W. H., and K~ong, A. (1994), "Covariance Structure of the Gibbs
Sampler with Applications to Comparisons of Estimators and Augmentation
Schemes," Biometrika, 81, 27-40.

[26] -(1995), "Covariance Structure and Convergence Rate of the Gibbs Sampler with
Various Scans," Journal of the Roriarl Statistical 8... .:. It; Series B, 57, 157-169.

[27] Liu, J. S. and Wu, Y. N. (1999), "Parameter Expansion for Data Augmentation,"
Journal of the American Statistical Association, 94, 1264-1274.

[28] Marchev, D. and Hobert, J. P. (2004), "Geometric ergodicity of van Dyk and Meng's
algorithm for the multivariate Student's t model," Journal of the American Statistical
Association, 99, 228-238.

[29] Mardia, K(., K~ent, J., and Bibby, J. (1979), M~ultivariate A,...el;;-.: London: Academic
press.

[:30] ?1n i-n, S. P. and Tweedie, R. L. (199:3), Iabrkov C'hain~s and Stochrastic St.,t..:1i;,
London: Springfer Verlagf.

[:31] -(1994), "Computable bounds for geometric convergence rates of Markov chains,"
The Annals of Applied P ~rol~.:l..7.;; 4, 981-1011.

[:32] Mira, A. and Geyer, C. J. (1999), "Ordering Monte Carlo Markov chains," Tech. Rep.
No. 6:32, School of Statistics, University of Minnesota.

[:33] Mykland, P., Tierney, L., and Yu, B. (1995), "Regeneration in Markov chain
Samplers," Journal of the American Statistical Asmsociation, 90, 2:3:341.

[:34] Nummelin, E. (1984) General Irreducible Iabrkov C'hain~s and Non-negative Op~era-
tors, London: Cambridge University Press.

[:35] Peskun, P. H. (197:3), "Optimum Monte Carlo sampling using Markov chains,"
Biometrikes, 60, 607-612.

[:36] R Development Core Team (2006), R: A Litr,:l;,rll. and Environment for Statistical
C'or,,ly-l.::l R Foundation for Statistical Computing, Vienna, Austria.

[:37] Retherford, J. R. (199:3) Hilbert Space: C'ompact Op~erators and the Trace theorem,
Cambridge University Press.

[:38] Robert, C. and Casella, G. (2004), M~onte Caerlo Statistical M~ethods, Springer, New
York.

[:39] Robert, C. P. (1995), "Simulation of truncated normal variables," Stratistic~s and

[40] Roberts, G. and Tweedie, R. (1999), "Bounds on regeneration times and convergence
rates for Markov chains," Stocheastic Processes and their Applications, 80, 221-229.
Corrigendum (2001) 91: 3:37-:338.

[41] Roberts, G. O. and Rosenthal, J. S. (1997), "Geometric ergodicity and hybrid Markov
chains," Electronic C'otmunications in Pr o~l~.:l..7.;; 2, 1:325.

[42] -(2001), :\1 I) 1:0v C'I !.1.4 and de-initializing processes," Scandinavian Journal of
Shetistic~s, 28, 489-504.

[4:3] -(2004), "General State Space Markov ('1. !.1.4 and MC' lC Algorithms," Pr o~l~.,t.7. /
Surveys, 1, 20-71.

[44] Roberts, G. O. and Tweedie, R. L. (2001), "Geometric L2 and L1 convergence are
equivalent for reversible Markov chains," Jourmal of Applied P l~r l~.:l..7.;; :38A, :37-41.

[45] Rosenthal, J. S. (1995), il!s...il. II!. .1, Conditions and Convergence Rates for Markov
C'I I!1, Monte Carlo," Journal of the American Statistical A~ssociation, 90, 558-566.

[46] Roy, V. and Hohert, J. P. (2007), "Convergence rates and .I-i nspinile; 1 standard errors
for MC1|LC algorithms for B li- -i Ils profit regression," Journal of the Roriarl Statistical
S .. .:. It; Series B, 69, 607-623.

[47] Rudin, W. (1991), Functional A,:tale;,.: McGraw-Hill, 2nd ed.

[48] Schervish, 31. J. and Carlin, B. P. (1992), "On the Convergence of Successive
Substitution CI .visplilrs Joureml of C'omputational and Grap~hical Stratistic~s, 1,
111-127.

[49] Tierney, L. (1994), 11 I) In~v chains for exploring posterior distributions (with
discussion)," The Annedls of Shetistic~s, 22, 1701-1762.

[50] -(1998), "A Note on Metropolis-Hastings K~ernels for General State Spaces," The
Annals of Applied Pr o~l~.:l..7.;; 8, 1-9.

[51] van Dyk, D. A. and Meng, X.-L. (2001), "The Art of Data Augmentation (with
Discussion) ," Jourmal of C'omputatiomel and Grap~hical Stratistic~s, 10, 1-50.

BIOGRAPHICAL SKETCH

Mr. Vivekananda Roy was born in 1980, in West Bengal, India. He spent his

childhood in his ancestral village, Bachhipur, where he went to school. After passing the

higher secondary examination, he moved to Calcutta in 1998. He received his bachelor's

degree in statistics from the University of Calcutta in 2001. He then joined the Indian

Statistical Institute, from where he received his Master of Statistics degree in 2003. He

joined the graduate program of the Department of Statistics at University of Florida in fall

2003. Upon graduation from UF, he will join the Department of Statistics at lowa State

University, as an Assistant Professor.

PAGE 1

1

PAGE 2

2

PAGE 3

3

PAGE 4

PAGE 5

page ACKNOWLEDGMENTS ................................. 4 LISTOFTABLES ..................................... 6 ABSTRACT ........................................ 7 CHAPTER 1INTRODUCTION .................................. 9 2MARKOVCHAINBACKGROUND ........................ 18 3BAYESIANPROBITREGRESSION ........................ 24 3.1Introduction ................................... 24 3.2GeometricConvergenceandCLTsfortheACAlgorithm .......... 27 3.3ComparingtheACandPX-DAAlgorithms ................. 35 3.4ConsistentEstimatorsofAsymptoticVariancesviaRegeneration ...... 37 4BAYESIANMULTIVARIATEREGRESSION 48 4.1Introduction ................................... 48 4.2ProofofPosteriorPropriety .......................... 50 4.3TheAlgorithms ................................. 58 4.3.1DataAugmentation ........................... 58 4.3.2HaarPX-DAAlgorithm ......................... 60 4.4GeometricErgodicityoftheAlgorithms .................... 61 5SPECTRALTHEOREMANDORDERINGOFMARKOVCHAINS ...... 71 5.1SpectralTheoryforNormalOperators .................... 71 5.2ApplicationofSpectralTheorytoMarkovChains .............. 81 APPENDIX:CHENANDSHAO'SCONDITIONS ................... 89 REFERENCES ....................................... 90 BIOGRAPHICALSKETCH ................................ 94 5

PAGE 6

Table page 3-1ResultsbasedonR=100regenerations ...................... 47 6

PAGE 7

7

PAGE 8

8

PAGE 9

PAGE 10

2exp 10

PAGE 11

2exp 2,theposteriordensitytakesthefollowingform(;jy)=1 2exp 2wherec2(y)isthemarginaldensityofygivenbyc2(y)=ZWZRdknYi=1Z10d 2exp 2dd;whereWRd(d+1) 2isthesetofddpositivedenitematrices.Inchapter4,weprovidenecessaryandsucientconditionsforc2(y)<1.Asinthepreviousexample,posteriorexpectationswithrespecttotheposteriordensity,(;jy),arenotavailableinclosed-form.Wenowdiscussdierentcomputationalmethodsthatcanbeusedtoapproximate( 1{1 ).Thesecomputationalmethodsarebroadlyoftwotypes,namely,numericalintegrationmethodsandsimulationbasedmethods.Ifthedimension,p,isnotlarge, 11

PAGE 12

PAGE 13

1{1 )bysimulatingaMarkovchainwithstationarydistribution.ThisisthebasicprincipleofMarkovchainMonteCarlo(MCMC)method.ThemostgeneralalgorithmforproducingMarkovchainswitharbitrarystationarydistributionistheMetropolis-Hastings(M-H)algorithm.AsimpleintroductiontotheM-HalgorithmisgiveninChibandGreenberg[ 7 ].Another 13

PAGE 14

4 ].Supposethep-dimensionalvectorin( 1{1 )canbewrittenas=(1;2;:::;p).ThesimplestGibbssampler(but,notthegeneralGibbssampler)requiresonetobeabletosimulatefromallunivariatefullconditionaldensitiesofi.e.,itisrequiredtosimulatefromtheconditionaldistributions,ijfj;j6=igfori=1;2;:::;p.ItisalsopossibletocreateahybridalgorithmwhichusesdierentversionsofM-HalgorithmtogetherwithGibbssamplertoconstructaMarkovchainwithstationarydistribution.Asourdiscussionsuggests,thereisaplethoraofMarkovchainswithstationarydistribution.InordertochoosebetweenMCMCalgorithms,weneedanorderingofMarkovchainshavingthesamestationarydistribution.InChapters2and5,wedescribedierentpartialorderingsofMarkovchains.LetfXjg1j=0denotetheMarkovchainassociatedwithanMCMCalgorithmthatisusedtoexplore.IffXjg1j=0isHarrisergodic(denedinChapter2),theergodictheoremimpliesthat,nomatterwhatthedistributionofthestartingvalue,X0, 14

PAGE 15

PAGE 16

AlbertandChib 's[ 1993 ]dataaugmentationalgorithmand LiuandWu 's[ 1999 ]PX-DAalgorithm.Westudytheconvergencerateofthesealgorithmsandprovetheexistenceofcentrallimittheorems(CLTs)forergodicaveragesunderasecondmomentcondition.WecomparethesetwoalgorithmsandshowthatthePX-DAalgorithmshouldalwaysbeusedsinceitismoreecientthantheotheralgorithminthesenseofhavingsmallerasymptoticvarianceinthecentrallimittheorem(CLT).Asimple,consistentestimatoroftheasymptoticvarianceintheCLTisconstructedusingregenerativesimulationmethods.InChapter4,weconsiderBayesianmultivariateregressionmodelswherethedistributionoftheerrorsisascalemixtureofnormals.Wenoticedbeforethatifthestandardnoninformativepriorisusedontheparameters(;),thenposteriorexpectationswithrespecttothecorrespondingposteriordensity,(;jy),arenotavailableinclosed-form.WedeveloptwoMCMCalgorithmsthatcanbeusedtoexplorethedensity(;jy).ThesealgorithmsarethedataaugmentationalgorithmandtheHaarPX-DAalgorithm.Wecomparethetwoalgorithmsandstudytheirconvergerates.Wealsoprovidenecessaryandsucientconditionsfortheproprietyoftheposteriordensity,(;jy).WhileinChapters3and4,weusedprobabilistictechniquestoanalyzedierentMCMCalgorithms,itispossibletotakeafunctionalanalyticapproachtostudyandcomparedierentMarkovchains.InChapter5,wegiveabriefoverviewofsomeresultsfromfunctionalanalysis.Inparticular,wediscussthespectraltheoremforbounded, 16

PAGE 17

17

PAGE 19

PAGE 20

PAGE 21

Rosenthal 's[ 1995 ]Theorem12showsthattheabovedriftandminorizationconditions,together,implythatisgeometricallyergodic.InChapter4,weemploydriftandminorizationconditionstoprovethegeometricergodicityofthedataaugmentationalgorithmusedinBayesianmultivariateStudent'stregressionproblem.Oneadvantageofprovinggeometricergodicityofbyestablishingtheabovedriftandminorizationconditionsisthatusing Rosenthal 's[ 1995 ]Theorem12,wealsocancalculateanupperboundofM(x)min( 2{1 ).Thisupperboundcanbeusedtocomputeanappropriateburn-inperiod(JonesandHobert[ 23 ],MarchevandHobert[ 28 ]).ThereareothermethodsofprovinggeometricergodicityofaMarkovchainthatdonotprovideanyquantitativeboundofM(x)min( 2{1 ).Wedescribeonesuchmethodnow.WewillassumethatXisequippedwithalocallycompact,separable,metrizabletopologywithB(X)astheBorel-eld.AfunctionV:X![0;1)issaidtobeunboundedocompactsetsifforevery>0,thelevelsetfx:V(x)giscompact.TheMarkovchainissaidtobeaFellerchainif,foranyopensetO2B(X),P(;O)isalower-semicontinuousfunction.Thefollowingpropositionisaspecialcaseof MeynandTweedie 's[ 1993 ]Lemma15.2.8.

PAGE 22

1 toestablishgeometricergodicityofMCMCalgorithmsusedinBayesianprobitregressionproblem.HobertandGeyer[ 15 ]employedProposition 1 toestablishthegeometricergodicityofGibbssamplersassociatedwithBayesianhierarchicalrandomeectsmodels.Noticethat,unlikeProposition 1 ,thedriftconditionin Rosenthal 's[ 1995 ]Theorem12doesnotrequirethedriftfunction,V,tobeunboundedocompactsets.Also, Rosenthal 's[ 1995 ]Theorem12doesnotneedtobeaFellerchain.ThedrivingforcebehindMCMCistheergodictheorem,whichissimplyaversionofthestronglawthatholdsforwell-behavedMarkovchains;e.g.,HarrisergodicMarkovchains.Indeed,supposethatf:X!RissuchthatRXjfjd<1anddeneEf=RXfd.Thentheergodictheoremsaysthattheaverage 41 ].(FormoreontheCLTinMCMC,seeChanandGeyer[ 5 ],MiraandGeyer[ 32 ],Jones[ 20 ]andJonesetal.[ 21 ].)ForathoroughdevelopmentofgeneralstatespaceMarkovchaintheory,seeNummelin[ 34 ]andMeynandTweedie[ 30 ].RobertsandRosenthal[ 43 ]providesaconcise,self-containeddescriptionongeneralstatespaceMarkovchains(alsoseeTierney[ 49 ]).AsmentionedinChapter1,foragivendistributionfunction,,therearelargenumberofMCMCalgorithmswithstationarydistribution.Onewaytoorderthese 22

PAGE 23

PAGE 24

6 ]providenecessaryandsucientconditionsonyandfxigni=1forc1(y)<1andtheseconditionsarestatedexplicitlyintheAppendix.Whentheseconditionshold,theposteriordensityofiswelldened(i.e.,proper)andisgivenby(jy)=1 AlbertandChib 's[ 1993 ]dataaugmentationalgorithm,whichwenowdescribe.LetXdenotethenpdesignmatrixwhoseithrowisxTiand,forz=(z1;:::;zn)T2Rn,let^=^(z)=(XTX)1XTz.Also,letTN(;2;w)denoteanormaldistributionwith 24

PAGE 25

(i) Drawz1;:::;znindependentlywithziTN(xTi;1;yi) (ii) Draw0Np^(z);(XTX)1AlbertandChib[ 1 ]hasbeenreferencedover350times,whichshowsthattheACalgorithmanditsvariantshavebeenwidelyappliedandstudied.ThePX-DAalgorithmofLiuandWu[ 27 ]isamodiedversionoftheACalgorithmthatalsosimulatesaMarkovchainwhoseinvariantdensityis(jy).AsingleiterationofthePX-DAalgorithmentailsthefollowingthreesteps: (i) Drawz1;:::;znindependentlywithziTN(xTi;1;yi) (ii) Drawg2Gamman 2Pni=1zixTi(XTX)1XTz2andsetz0=(gz1;:::;gzn)T Draw0Np^(z0);(XTX)1NotethattherstandthirdstepsofthePX-DAalgorithmarethesameasthetwostepsoftheACalgorithmso,nomatterwhatthedimensionof,thedierencebetweentheACandPX-DAalgorithmsisjustasingledrawfromtheunivariategammadistribution.Fortypicalvaluesofnandp,theeortrequiredtomakethisextraunivariatedrawisinsignicantrelativetothetotalamountofcomputationneededtoperformoneiterationoftheACalgorithm.Thus,thetwoalgorithmsarebasicallyequivalentfromacomputationalstandpoint.However,LiuandWu[ 27 ]andvanDykandMeng[ 51 ]bothprovideconsiderableempiricalevidencethatautocorrelationsdiedownmuchfasterunderPX-DAthanunderAC,whichsuggeststhatthePX-DAalgorithm\mixesfaster"thantheACalgorithm.(LiuandWu[ 27 ]alsoestablishedatheoreticalresultalongtheselines-seetheproofofourCorollary 1 .) 25

PAGE 26

3{1 )isbyestablishinggeometricergodicityoffjg1j=0.Inthischapter,weprovethattheMarkovchainsunderlyingtheACandPX-DAalgorithmsbothconvergeatageometricratewhichimpliesthattheCLTin( 3{1 )holdsforeveryf2L2(jy);thatis,foreveryfsuchthatRRpf2()(jy)d<1.WealsoestablishthatPX-DAistheoreticallymoreecientthanACinthesensethattheasymptoticvarianceintheCLTunderthePX-DAalgorithmisnolargerthanthatundertheACalgorithm.Regenerativemethodsareusedtoconstructasimple,consistentestimatoroftheasymptoticvarianceintheCLT.Asanillustration,weapplyourresultsto vanDykandMeng 's[ 2001 ]lupusdata.Inthisparticularexample,theestimatedasymptoticrelativeeciencyofthePX-DAalgorithmwithrespecttotheACalgorithmisabout65.Hence,eventhoughtheACandPX-DAalgorithmsareessentiallyequivalentintermsofcomputationalcomplexity,hugegainsineciencyarepossiblebyusingPX-DA. 26

PAGE 27

3.2 and 3.3 ,respectively.InSection 3.4 wederiveresultsthatallowfortheconsistentestimationofasymptoticvariancesviaregenerativesimulation.

PAGE 28

34 ,Theorem3.8].Supposehisabounded,harmonicfunction.SincetheACalgorithmis-irreducibleandhasaninvariantprobability 28

PAGE 29

34 ,Proposition3.13].Thus,thereexistsasetNwith(N)=0suchthath()=cforall2 2 itfollowsthatergodictheoremholdsforit.Thefollowingtheoremisthemainresultofthissection. Proof. 1 .WehaveshownthatACalgorithmis-irreducibleandaperiodic,wheredenotetheLebesguemeasureonRp.So,ifisamaximalirreducibilitymeasurefortheMarkovchainunderlyingtheACalgorithm,then.Conversely,if(A)=0,thenKm(;A)=0forall2Rpandallm2N,whichimpliesthat(A)=0anditfollowsthat.Hence,.Sincethesupportofobviouslyhasnon-emptyinterior,itfollowsthatthesupportofamaximalirreducibilitymeasurefortheACalgorithmhasnon-emptyinterior.WenowdemonstratethattheMarkovchainassociatedwiththeACalgorithmisaFellerchain.LetandOdenoteapointandanopensetinRp,respectively.Assumethatfmg1m=1;isa(deterministic)sequenceinRpwithm6=suchthatm!asm!1.TwoapplicationsofFatou'sLemmainconjunctionwiththefactthat(zj;y)iscontinuousinyieldliminfm!1K(m;O)ZOliminfm!1k(0jm)d0=ZOliminfm!1"ZRm(0jz;y)(zjm;y)dz#d0

PAGE 30

1 withdriftfunctionV()=(X)T(X).RecallthatXisassumedtohavefullcolumnrank,p,andhenceXTXispositivedenite.Thus,foreach>0,thesetf2Rp:V()g=f2Rp:T(XTX)giscompactsothefunctionVisunboundedocompactsets.Now,notethat(KV)()=ZRpV(0)k(0j)d0=ZRp"ZRnV(0)(0jz;y)d0#(zj;y)dz=ZRpEV(0)jz;y(zj;y)dz=EnEV(0)z;y;yo;where,asthenotationsuggests,theexpectationsinthelasttwolinesarewithrespecttotheconditionaldensities(0jz;y)and(zj;y).Recallthat(0jz;y)isap-dimensionalnormaldensityand(zj;y)isaproductoftruncatednormals.Evaluatingtheinsideexpectation,wehaveEV(0)z;y=E(0)TXTX0z;y=tr(XTX(XTX)1)+zTX(XTX)1(XTX)(XTX)1XTz=p+zTX(XTX)1XTzp+zTz;

PAGE 31

19 ]implythatifUTN(;1;1)then,E(U2)=1+2+() ();where()withonlyasingleargumentdenotesthestandardnormaldensityfunction;thatis,(v)isequivalentto(v;0;1).Similarly,ifUTN(;1;0)then,E(U2)=1+2() 1():ItfollowsthatEz2ij;y=8><>:1+(xTi)2+(xTi)(xTi) (xTi)ifyi=11+(xTi)2(xTi)(xTi) 1(xTi)ifyi=0:Amorecompactwayofexpressingthisisasfollows: Ez2ij;y=1+(wTi)2(wTi)(wTi) 1(wTi);(3{2)wherewiisdenedintheAppendix.Hence,wehave (KV)()=EnEV(0)z;y;yop+n+nXi=1(wTi)2nXi=1(wTi)(wTi) 1(wTi):(3{3)Recallthatthegoalistoshowthat(KV)()V()+Lforall2Rp.Itfollowsfrom( 3{3 )that(KV)(0)p+n.Wenowconcentrateon2Rpnf0g.WebeginbyconstructingapartitionofthesetRpnf0gusingthenhyperplanesdenedbywTi=0.Forapositiveintegerm,deneNm=f1;2;:::;mg.LetA1;A2;:::;A2ndenoteallthesubsetsofNn,and,foreachj2N2n,deneacorresponding 31

PAGE 32

[2nj=1Sj=Rpnf0g,and 5 areinforce,thereexiststrictlypositiveconstantsfaigni=1suchthata1wT1+a2wT2++anwTn=0:Therefore, 3{4 )impliesthattheremustalsoexistani06=isuchthatwTi0andwTihaveoppositesigns.Thus,AjandAjarebothnonempty.NowdeneC=j2N2n:Sj6=;.Foreachj2C,deneRj()=Pi2Aj(wTi)2 11 ,p.175].Also,itisclearthatif 32

PAGE 33

1(u);thenM2(0;1).Fixj2C.Itfollowsfrom( 3{3 )andtheresultsconcerningMillsratiothatforall2Sj,wehave(KV)()p+n+nXi=1(wTi)2Xi2Aj(wTi)(wTi) 1(wTi)Xi2 1(wTi)p+n+nXi=1(wTi)2+Xi2Aj(wTi)(wTi) 1(wTi)Xi2 1(wTi)p+n+nXi=1(wTi)2+nMXi2

PAGE 35

PAGE 36

10 ]forbackgroundonleftgroupactionsandmultipliers.)LetZdenotethesubsetofRninwhichzlives;i.e.,ZistheCartesianproductofnhalf-lines(R+andR),wheretheithcomponentisR+ifyi=1andRifyi=0.Fixz2Z.ItiseasytoseethatStep2ofthePX-DAalgorithmisequivalenttothetransitionz!gzwheregisdrawnfromadistributiononGhavingdensityfunction 17 ].First,theirProposition3showsthatR(z;dz0)isreversiblewithrespectto(zjy)anditfollowsthatk(0j)isreversiblewithrespectto(jy).WenowusethefactthattheACalgorithmisgeometricallyergodictoestablishthatthePX-DAalgorithmenjoysthispropertyaswell. Proof. 25 32 ].DenotethenormsoftheseoperatorsbykKkandkKk.Ingeneral,areversible,HarrisergodicMarkovchainisgeometricallyergodicifandonlyifthenormoftheassociatedMarkovoperator 36

PAGE 37

41 44 ].ByTheorem 1 ,theACalgorithmisgeometricallyergodicandconsequentlykKk<1.ButLiuandWu[ 27 ]showthatkKkkKk[seealso 17 ,Theorem4]andhencekKk<1,whichimpliesthatthePX-DAalgorithmisalsogeometricallyergodic. WehavenowshownthattheMarkovchainsunderlyingtheACandPX-DAalgorithmsarebothreversibleandgeometricallyergodicandhencebothhaveCLTsforallf2L2(jy).WenowuseanotherresultfromHobertandMarchev[ 17 ]toshowthatthePX-DAalgorithmisatleastasecientastheACalgorithm. Proof. HobertandMarchev 's[ 2008 ]Theorem4. Inordertoputourtheoreticalresultstouseinpracticetocomputevalidasymptoticstandarderrors,werequireaconsistentestimatoroftheasymptoticvarianceandthisisthesubjectofthenextsection. 3 .Ofcourse,themarginalchainfjg1j=0hastheMarkovtransitiondensityk(0j)nomatterwhichversionoftheMarkovtransitiondensitywechooseforthejointchain.Whilethisisobviousforthechaincorresponding~~k,itcanbeeasilyshownfor~kbyconsideringtwoconsecutivesteps 37

PAGE 38

42 ]canbeusedtoshowthatthejointchainfj;zjg1j=0inheritsgeometricergodicityfromitsmarginalchainfjg1j=0.Notethatfj;zjg1j=0isthechainthatisactuallysimulatedwhentheACalgorithmisrun(wejustignorethezjs).Supposewecanndafunctions:RpRn![0;1],whoseexpectationwithrespectto(;zjy)isstrictlypositive,andaprobabilitydensityd(0;z0)onRpRnsuchthatforall(0;z0);(;z)2RpRn,wehave ~k0;z0j;zs(;z)d(0;z0):(3{6)Thisiscalledaminorizationcondition[ 22 30 43 ]anditcanbeusedtointroduceregenerationsintotheMarkovchaindrivenby~k.Theseregenerationsarethekeytoconstructingasimple,consistentestimatorofthevarianceintheCLT.Afterexplainingexactlyhowthisisdone,wewillidentifysanddforbothACandPX-DA.Equation( 3{6 )allowsustorewrite~kasthefollowingtwo-componentmixturedensity ~k0;z0j;z=s(;z)d(0;z0)+1s(;z)r0;z0j;z;(3{7) 38

PAGE 39

1s(;z);whens(;z)<1(anddenedarbitrarilywhens(;z)=1).InsteadofsimulatingtheMarkovchainfj;zjg1j=0intheusualwaythatalternatesbetweendrawsfrom(zj;y)and(jz;y),wecouldsimulatethechainusingthemixturerepresentation( 3{7 )asfollows.Supposethecurrentstateis(j;zj)=(;z).First,wedrawjBernoullis(;z).Thenifj=1,wedraw(j+1;zj+1)fromd,andifj=0,wedraw(j+1;zj+1)fromtheresidualdensity.The(random)timesatwhichj=1correspondtoregenerationsinthesensethattheprocessprobabilisticallyrestartsitselfatthenextiteration.Morespecically,supposewestartbydrawing(0;z0)d.Theneverytimej=1,wehave(j+1;zj+1)dsotheprocessis,ineect,startingoveragain.Furthermore,the\tours"takenbythechaininbetweentheseembeddedregenerationtimesareiid,whichmeansthatstandardiidtheorycanbeusedtoanalyzetheasymptoticbehaviorofergodicaverages,therebycircumventingthedicultiesassociatedwithanalyzingaveragesofdependentrandomvariables.Formoredetailsandsimpleexamples,seeMyklandetal.[ 33 ]andHobertetal.[ 16 ].Inpractice,wecanevenavoidhavingtodrawfromr(whichcanbeproblematic)simplybydoingthingsinaslightlydierentorder.Indeed,giventhecurrentstate(j;zj)=(;z),wedraw(j+1;zj+1)intheusualway(thatis,bydrawingfrom(zj;y)and(jz;y))afterwhichwe\llin"avalueforjbydrawingfromtheconditionaldistributionofjgiven(j;zj)and(j+1;zj+1),whichisjustaBernoullidistributionwithsuccessprobabilitygivenby ~k(j+1;zj+1jj;zj):(3{8) 39

PAGE 40

N=1 16 ]areapplicableandimplythat,aslongasthereexistsan>0suchthatEjf()j2+jy<1,then N2: 3{9 )isslightlydierentfromtheCLTdiscussedearlier,whichtakestheformp

PAGE 41

PAGE 42

(2)p 2exp1 2^(z)TXTX^(z):where^(z)=(XTX)1XTzThus,s(z)="inf2D(jz;y) 2^(z)TXTX^(z)2TXTX^(z) 2^(z)TXTX^(z)2TXTX^(z)="exp1 2zTX(XTX)1XTz 2zTX(XTX)1XTzinf2Dexp(zz)TX="exp1 2zTX(XTX)1XTz 2zTX(XTX)1XTzexppXi=1citiIR+(ti)+ditiIR(ti)

PAGE 43

PAGE 44

PAGE 45

~k(j+1;(zj+1;gj+1)jj;(zj;gj))=inf2D(jgjzj;y) 2g2jzTjX(XTX)1XTzj 2zTX(XTX)1XTzexppXi=1cit(j)iIR+(t(j)i)+dit(j)iIR(t(j)i)exp1 2zTX(XTX)1XTz 2g2jzTjX(XTX)1XTzjexp(gjzjz)TXj+1ID(j+1)=exp(pXi=1cit(j)iIR+(t(j)i)+dit(j)iIR(t(j)i)t(j)ij+1;i)ID(j+1)wheret(j)T=(gjzjz)TX.Theorem 2 statesthattheasymptoticvarianceintheCLTforthePX-DAalgorithmisnolargerthanthatfortheACalgorithm;i.e,0<2f;k2f;k<1forallf2L2(jy).However,weknowfromRemark 1 thattheregenerativemethodisbasedonaslightlydierentCLTwhoseasymptoticvariancehasanextrafactorinvolvingthesmallfunction,namelyE(s()jy),fromtheminorizationcondition.AlthoughthesmallfunctionsinthetwominorizationconditionsthatwederivedfortheACandPX-DAalgorithmsareslightlydierent,E(s(:)jy)remainsexactlythesameasshown 45

PAGE 46

2(n2)=2(n=2)ZRnZR+z0T(IH)z0n vanDykandMeng 's[ 2001 ]lupusdata,whichconsistsoftriples(yi;xi1;xi2),i=1;:::;55,wherexi1andxi2arecovariatesindicatingthelevelsofcertainantibodiesintheithindividualandyiisanindicatorforlatentmembranouslupusnephritis(1forpresenceand0forabsence).vanDykandMeng[ 51 ]consideredthemodelPr(Yi=1)=0+1xi1+2xi2;withaatprioron.Weusedalinearprogram(thatisdescribedintheAppendix)toverifythat ChenandShao 's[ 2000 ]necessaryandsucientconditionsforproprietyaresatisedinthiscase. 46

PAGE 47

39 ]togenerateone-sidedtruncatednormalrandomvariables.WeranACandPX-DAforR=100regenerationseach.Thistook1,223,576iterationsforACand1,256,677iterationsforPX-DA.WeusedthesimulationstoestimatetheposteriorexpectationsoftheregressionparametersandtheresultsareshowninTable 3-1 .(ResultsinChenandShao[ 6 ]implythatthereexists>0suchthatEjjj2+jy<1forj2f0;1;2g.)ItisstrikingthattheestimatedasymptoticvariancesfortheACalgorithmareallatleast65timesaslargeasthecorrespondingvaluesforthePX-DAalgorithm.Theseestimatessuggestthat,inthisparticularexample,theACalgorithmrequiresabout65timesasmanyiterationsasthePX-DAalgorithmtoachievethesamelevelofprecision.(Weactuallyrepeatedtheentireexperimentseventimesandtheestimatesof2f;k=2f;krangedbetween40and145.) Table3-1. ResultsbasedonR=100regenerations ACAlgorithmPX-DAAlgorithm Parameter estimates.e.^f;k=p -3.0180.012 66.61 6.9160.023 66.92 3.9820.015 63.1 47

PAGE 48

2exp 4{1 )asy=X+"wherey=(y1;:::;yn)Tisthendmatrixofobservations,X=(x1;x2;:::;xn)Tisthenkmatrixofcovariatesand"=("1;:::;"n)Tisthendmatrixoferrorvariables.Thelikelihoodfunctionfortheregressionmodelin( 4{1 )isgivenbyf(yj;)=nYi=1Z10d 2exp

PAGE 49

2.Theposteriordensitytakesthefollowingform 2isthesetofddpositivedenitematrices.FernandezandSteel[ 12 ]provedthatc2(y)<1ifandonlyifnd+k.Insection 4.2 ,wegiveanalternativeproofoftheposteriorpropriety.Abyproductofourproofisamethodofexactsamplingfrom(;jy)intheparticularcasewhennisexactlyd+k.Throughoutthischapterweassumethatnd+k.Wealsoassumethatthecovariatematrix,X,isoffullcolumnranki.e.,r(X)=k.Theposteriordensityin( 4{3 )isintractableinthesensethatposteriorexpectationsarenotavailableinclosed-form.Also,ourexperienceshowsthatitisdiculttodevelopausefulprocedureformakingi.i.d.drawsfrom(;jy).InthischapterwefocusonMCMCmethodsforexploringtheposteriordensityin( 4{3 ).Wedevelopadataaugmentation(DA)algorithmfor(;jy)insection 4.3.1 .IthasbeennoticedintheliteraturethatthestandardDAalgorithmoftensuersfromslowconvergence[ 51 ].EmpiricalandtheoreticalstudieshaveshownthatalternativealgorithmsthataremodiedversionsofthestandardDAalgorithmsuchastheHaarPX-DAalgorithmandthemarginalaugmentationalgorithmoftenprovidehugeimprovementoverthestandardDAalgorithm(LiuandWu[ 27 ],vanDykandMeng[ 51 ],RoyandHobert[ 46 ],HobertandMarchev[ 17 ]).Insection 4.3.2 ,wedeveloptheHaarPX-DAalgorithmfortheposteriordensityin( 4{3 ).Wethenspecializetothecasewhentheerrors,"i's,haveaStudent'stdistributioni.e.,themixingdistribution,H(),in( 4{2 )isaGamma( 49

PAGE 50

17 ]wealsoconcludethattheHaarPX-DAalgorithmisatleastasecientasthedataaugmentationalgorithminthesensethattheasymptoticvariancesinthecentraltheoremunderHaarPX-DAalgorithmareneverlargerthanthoseundertheDAalgorithm.SomeoftheseresultsaregeneralizationsofresultsfromvanDykandMeng[ 51 ]andMarchevandHobert[ 28 ]whoconsideredthespecialcasewheretherearenocovariatesintheregressionmodel( 4{1 ),i.e.,X=(1;1;:::;1)TandHisGamma( 4.3 ,wedescribetheDAandtheHaarPX-DAalgorithms.Inthelastsection,wecomparethetwoalgorithmsandprovethatboththealgorithmsconvergeatageometricrate. 50

PAGE 51

4{1 )and( 4{2 )itfollowsthat (2)nd 2:

PAGE 52

2 ,Chapter17] 2 ]Supposeisanmrmatrix,andaremmandrrnon-negativedenitematrices.WesaythatZhasamatrixnormaldistributionwithparameters,andifZisanmrrandommatrixhavingmomentgeneratingfunctiongivenbyMZ(t)=exptr(Tt)+1 2tr(tTt)andwewriteZNm;r(;;).InthiscasewehaveE(Z)=.Moreover,ifandarepositivedenitematricesthenZhasthefollowingdensityfunctionfZ(z)=1 (2)mr=2jjr=2jjm=2exp1 2tr1(z)1(z)T:Sincer(X)=k,itfollowsthatXTQ1Xandhenceisap.d.matrix.Thus, 2tr1yTQ1yT1)jQjd 2

PAGE 53

2 2exp(1 2tr1yTQ1yT1)jjd 2IQ1 2X(XTQ1X)1XTQ1 2Q1 2y andsinceIQ1 2X(XTQ1X)1XTQ1 2isanidempotentmatrix,itimpliesthatyTQ1yT1isapositivesemi-denitematrix.NowweprovethatyTQ1yT1isap.d.matrixbyshowingthatjyTQ1yT1j6=0(withprobabilityone).Letbethen(d+k)augmentedmatrix(X:y).Then,TQ1=264XTyT375Q1Xy=264XTQ1XXTQ1yyTQ1XyTQ1y375:Therefore, Sincer(X)=k,weknowthatXTQ1Xisap.d.matrixandhencejj>0.Alsosincend+k,thed+kcolumnsofarelinearlyindependentwithprobabilityonebecausetheprobabilityofndimensionalrandomcolumnvectorsofylyinginanylinearsubspaceofRnwithdimensionn1iszero(withrespecttotheLebesguemeasureonRdn). 53

PAGE 54

4{8 )itfollowsthatjyTQ1yT1j>0.Tointegratetheexpressionin( 4{6 )withrespectto,weusethefollowingdenitionofInverseWishartdistribution. 29 ,p.85]Letbeappp.d.matrix.ThenforsomempthepprandommatrixWissaidtohaveaInverseWishartdistributionwithparametersmandifthep.d.fofW(withrespecttoLebesguemeasureonRp(p+1) 2,restrictedtothesetwhereW>0)isgivenbyf(W;m;)=jWjm+p+1 2exptr(1W1) 2 4jjm 2(m+1i))andwewriteWIWp(m;).Henceifnd+ki.e.,nkdbytheabovedenitionofInverseWishartdistribution,wehave =2d(nk) 2d(d1) 4Qdi=1(1 2(nk+1i)) (2)d(nk) 2yTQ1yT1nk 2(nk+1i)) 4yTQ1yT1nk 4{8 )weget 4jTQ1jd 4jjdnYj=1h(qj): 54

PAGE 55

4jjd;which,ofcourse,isanitenumber.So,wehaveprovedthatintheparticularcasewhenn=d+k,theposteriordistributionisproperwithprobabilityone.Then,anapplicationofLemma2ofMarchevandHobert[ 28 ]showsthatfor-almostallytheposteriorisproperfornd+k.Wewillnowshowthattheposteriordistributionisimproperwhenn
PAGE 56

2:uii>0;8i=1;:::;dg:NoticethatthecolumnsofLformabasisofRd.Sothereexistsconstantsb1;b2;:::;bdwithbi6=0forsomeisuchthatx0=Pdi=1bili.Supposei0=minfi2f1;2;:::;dg:bi6=0g.NowconsiderthetransformationL=(l1;l2;:::;ld)!O=(o1;o2;:::;od)whereoi=li8i6=i0andoi0=x0i.e.,O=LA,whereA=0BBBBBBBBBBBBBB@10000100...00bi00...00bd11CCCCCCCCCCCCCCA:ThentheJacobianofthetransformationisjjAjdj=jbi0jd[ 29 ,p.36].Notethatli0=1 2dXi=1lTilidYi=1lnkiiidL=jbi0jd 2Xi6=i0oTioi+1 2Xii01+b2i 2:vii>0;8i6=i0andbi0vi0i0>0g:ByFubini'stheoremwecanrearrangetheorderofintegration.Noticethatoi0doesnotappearintheexponentialtermandtheonlyterminvolvingoi0i0intheaboveintegralis 56

PAGE 57

AswementionedintheintroductionthatFernandezandSteel[ 12 ]gaveaproofofproprietyoftheposteriordensity(;jy).Abyproductofouralternativeproofisamethodofexactsamplingfrom(;jy)intheparticularcasewhennisd+k.Wedescribethemethodnow.Let(q;;jy)bethejointposteriordensityof(q;;)givenby(q;;jy)=f(y;qj;)(;) 4{5 )itiseasytoseethat 4{10 )weknowthat(qjy)=Qni=1h(qi).So,inthiscaseanexactdrawfrom(;jy)canbemadeusingthefollowingthreesteps: (i) Drawq1;q2;:::qnindependentlywhereqih(qi). (ii) DrawIWd[nk;(yTQ1yyTQ1X(XTQ1X)1XTQ1y)1] (iii) DrawTNd;k(yTQ1X(XTQ1X)1;;(XTQ1X)1) 57

PAGE 58

36 ])havefunctionsforgeneratingrandommatricesfromtheInverseWishartdistribution.OnewaytogenerateZNm;r(;;)istorstindependentlydrawZiindNr(i;)whereTiistheithrowoffori=1;:::;m.ThentakeZ=1 2266666664ZT1ZT2...ZTm377777775: 4{3 ).WerstdeveloptheDAalgorithminSection 4.3.1 usingthelatentdataq=(q1;q2;;qn)andthejointposteriordensityof(q;;jy).WethenderivetheHaarPX-DAalgorithminsection 4.3.2 .Inthespecialcase,whentheobservations,yi's,areassumedtobefromamultivariateStudent'stlocation-scalemodel,vanDykandMeng[ 51 ]developedthemarginalaugmentationalgorithm,whichisamodiedversionofthestandarddataaugmentationalgorithm,forthedensity( 4{3 ).HobertandMarchev[ 17 ]haveshownthatwhenthegroupstructureexploitedbyLiuandWu[ 27 ]exists,marginalaugmentationalgorithm(withleft-Haarmeasurefortheworkingprior)isexactlysameasthe LiuandWu 's[ 1999 ]HaarPX-DAalgorithm.Insection 4.3.2 weshowthatasimilargroupstructurecanbeestablishedforanalyzingtheposteriordensityin( 4{3 )andsomarginalaugmentationalgorithmisthesameastheHaarPX-DAalgorithminourcase.

PAGE 59

4.2 ,weknowthatconditionallyj;q;yfollowsaMatrixNormaldistributionandtheconditionaldistributionofjq;yisanInverseWishartdistribution.Conditionalon(;;y)qi'sareindependentwithqij;;yindh(qi)qd 2:

PAGE 60

HobertandMarchev 's[ 2008 ]Section4.3.FromSection 4.2 ,weknowthat(qjy)/jXTQ1Xjd q:d q:nk g=c3Z10nYi=11 :Inordertoshowthatm(q)<1,supposex1 ;x>0andconsiderthestandardnoninformativeprior,1=,forthescaleparameter,.ThenthecorrespondingposteriordistributionispropersinceZ101 d =1 )=1

PAGE 61

28 ]itfollowsthatR10Qni=11 <1forallx1;x2;:::;xn>0andn1sinceQni=11 .Hence,itfollowsthatm(q)<1.ConsiderthefollowingunivariatedensityonR+eq(g)=gn1 (i) Drawq(qj0;0;y) (ii) Drawgeq(g)andsetq0=gq Draw;(;jq0;y)TheMarkovtransitiondensityofthePX-DAalgorithmcanbewrittenas~k(;j0;0)=ZRn+ZRn+(;jq0;y)R(q;dq0)(qj0;0;y)dq;whereR(q;dq0)istheMarkovtransitiondensityinducedbythestep2,thattakesq!q0=gq.InthespecialcasewhenwehavemultivariateStudent'stdata,itiseasytoseethatthedensityeq(g)isGamma(n 61

PAGE 62

jPj:SincePisap.d.matrix,jPj>0.Similarly,sincePxxTisap.s.d.matrix,jPxxTj0.ThenfromtheaboveidentityitfollowsthatxTP1x1. ThefollowinglemmaestablishesthedriftinequalityfortheDAalgorithm. +d2V(0;0)+n+dk +d2:

PAGE 63

Tocalculatetheaboveconditionalexpectationsweneedthecorrespondingconditionaldistributions.Therequiredconditionaldistributions,asderivedintheprevioussections,arethefollowingTj;q;yNd;k(T;;)jq;yIWdnk;yTQ1yT11andqij;;yind+d 2i=1;2;:::;n:Startingwiththeinnermostexpectation,wehaveE(V(;)j;q;y)=EnXi=1yTi1yi2nXi=1yTi1Txi+nXi=1xTi1Txi;q;y=nXi=1yTi1yi2nXi=1EyTi1Txi;q;y+nXi=1ExTi1Txi;q;y:Tocalculatetheaboveexpectationsweusethefollowingpropertyofmatrixnormaldistribution.LetZNm;r(;;).IfCandDaretwomatricesofdimensionnmandrsrespectivelythenCZDNn;s(CD;CC0;D0D)[ 2 ,chap17].So,1 2T;q;yNd;k(1 2T;I;).Notethat,1T=1 2TT(1 2T.Therefore,1T;q;yWk(d;;1T),wherebyVWp(m;;)wemeanthatVhaspdimensionalnoncentralWishartdistributionwithmdegreesoffreedom,covariancematrixandwithnoncentralityparameter[ 2 ,chap17].Inthiscase,E(V)=m+.If=0,wesaythatVhasacentralWishartdistributionandwewriteVWp(m;).So,wehave 63

PAGE 64

4{12 ),weusethefactthatifXIWp(m;),thenX1Wp(m;)andE(X1)=m[ 29 ,p.85].Thus,E[E(V(;)j;q;y)jq;y]=nXi=1(yiTxi)TE1jq;y(yiTxi)+dnXi=1xTixi=(nk)nXi=1(yiTxi)T(yTQ1yT1)1(yiTxi)+dnXi=1xTixi:

PAGE 65

1 yields (yiTxi)T(yTQ1yT1)1(yiTxi)1 1 gives Therefore,wegetE(V(;)jq;y)(n+dk)nXi=11 2i=1;2;:::;n:So,usingthefactthatifwGamma(a;b)thenE(1 a1,nally,wehaveE(V(;)j0;0)n+dk +d2nXi=1(yi0Txi)T(0)1(yi0Txi)+;whichprovesthelemma. Thefollowinglemmaestablishesanassociatedminorizationcondition. 65

PAGE 66

llog1+l and(a;b;x)denotestheGamma(a,b)densityevaluatedatthepointx. Proof. 2nYi=1inf(0;0)2S+d 2nYi=1inf(0;0)2Si+d 2:Thenby Hobert 's[ 2001 ]Lemma1,itstraightforwardlyfollowsthat(qj0;0;y)nYi=1g(qi)8(0;0)2S:

PAGE 67

+d2isstrictlylessthan1.Westatethisinthefollowingtheorem. +d2<1i.e.,n<+k2. 's[ 2008 ]Proposition6showsthattheHaarPX-DAalgorithmisatleastasecient(ineciencyordering)astheDAalgorithm.UsingsimilarargumentsasinCorollary1,wecanshowthatgeometricergodicityoftheDAalgorithmimpliesthatoftheHaarPX-DAalgorithm.Hencewehavethefollowingcorollary. 2 matcheswith MarchevandHobert 's[ 2004 ]Theorem10inthecasewhenk=1.WealsocanprovethegeometricergodicityoftheHaarPX-DAalgorithmbydirectlyestablishingadriftandminorizationconditionforit.WeactuallycanusethesamedriftfunctionV(;)=Pni=1(yiTxi)T1(yiTxi)toestablishadriftconditionfortheHaarPX-DAalgorithm. (n2)(+d2)V(0;0)+n(n+dk) (n2)(+d2):

PAGE 68

4.3.2 weknowthatq0=gqwheregjqGamma(n 4{13 )andthenstraightforwardalgebrashowsthatE(V(;)jq0;y)=(nk) gnXi=1xTixi:So,E(V(;)jq;y)=(nk)nXi=1(yiTxi)T(yTQ1yT1)1(yiTxi)+dnXi=1xTixiE1 4{14 )and( 4{15 )itthenfollowsthatE(V(;)jq;y)(n+dk) 2i=1;2;:::;n;wehaveE(V(;)j0;0)(n+dk) +d2nXi=1Xj6=i+(0Txiyi)T(0)1(0Txiyi)

PAGE 69

OurminorizationconditionfortheDAalgorithmstraightforwardlygeneralizestoaminorizationconditionforthePX-DAalgorithm. 3 .Together,Lemmas 4 and 5 provethefollowingtheorem. (n2)(+d2)<1.AsacorollaryofTheorem 4 (seeCorollary 2 )weknowthatthePX-DAalgorithmisgeometricallyergodicifn+dk +d2<1.AtrstitmightappearthatTheorem 5 isabetterresultthanCorollary 2 .But,wenowshowthatitcanneverhappenthatbothofthefollowinginequalitiesholdtogether +d21;(n+dk)(+d)(n1) (n2)(+d2)<1:(4{16)Notethat(n+dk)(+d)(n1) (n2)(+d2)
PAGE 70

4{17 )holdsifandonlyifn+k2andn<1+2 4 andTheorem 5 areupshotsoftheparticulardriftfunctionV(;)andtheinequalities( 4{14 )and( 4{15 )thatweusedtoprovethedriftconditions.Inouropinion,tomakesubstantialimprovementofTheorem 4 andTheorem 5 ,weeitherhavetoconsideradierentdriftfunctionorresorttoaaltogetherdierenttechniqueofprovinggeometricergodicityotherthanestablishingdriftandminorizationcondition.Eitherofthesetwowouldrequireustostartfromscratch. 70

PAGE 71

(h;g) 2.Iftheresultingnormedspaceiscomplete,itiscalledaHilbertspace.LetHbeaHilbertspaceoverCandT:H!Hbealineartransformation,calledalinearoperator.TheoperatorTissaidtobeboundedifthereexistsanM>0suchthatkThkMkhkforallh2H,wherethenorm,kk,isasdenedabove.LetB(H)bethecollectionofalllinear,boundedoperatorsfromHintoH.ForT2B(H),denetheoperatornormofTbykTk=supkThk khk:h2H;h6=0:ItiseasytoseethatkTk=supfkThk:h2H;khk1g=supfkThk:h2H;khk=1g:

PAGE 72

47 ,p.311].Somepropertiesofadjointoperatorsarelistedbelow:kTk=kTkT=TkTTk=kTk2andifSalsobelongstoB(H)then(ST)=TS:Wenowcandenedierenttypesofoperators.AnoperatorT2B(H)issaidtobe 4 weknowthat(h;Th)= (Th;h). 72

PAGE 73

47 ,p.310]. Proof. SupposeT2B(H).TheoperatorTisinvertibleifthereisanS2B(H)suchthatTS=I=ST.TheoperatorSiscalledtheinverseofTandwewriteS=T1.SupposeR(T)denotestherangeofT,i.e.,R(T)=fTh:h2Hg.ThenTisinvertibleifandonlyifR(T)=HandTisone-to-one.Thespectrum,(T),oftheoperatorT2B(H)isdenedasfollows(T)=f:TIisnotinvertibleg:Thus2(T)ifandonlyifatleastoneofthefollowingtwostatementsistrue: (i) TherangeofTIisnotallofHi.e.,TIisnotonto. (ii) TheoperatorTIisnotone-to-one. 73

PAGE 74

9 ,p.83]. n: Proof.

PAGE 75

4 tellsusthat(T)isaboundedsetforanyboundedlinearoperatorTbecause(T)f2C:jjkTkg: Proof. ItcanalsobeprovedthatforT2B(H);H6=f0g,thespectrumofT,(T),isnotempty[ 47 ,p.253].Thespectralradius,rT,ofTisdenedasrT=supfjj:2(T)g:

PAGE 76

PAGE 77

PAGE 78

PAGE 79

13 isofcourseaboundedoperatoranditcanbeprovedthatkf(T)ksupfjf()j:2(T)g.Wealsohavethefollowingtheoremwhichdeterminesthespectrumoftheoperatorf(T). 37 ]givesanice,concisedescriptionofcompactoperators.ForacompactoperatorT,thespectrum,(T),isatmostcountable,andhasatmostonelimitpoint,namely,0.Also,anynon-zerovalueinthespectrumisnecessarilyaneigenvalueofT.Westatetheseresultsinthefollowingtheorem([ 8 ]p.214,[ 9 ]chapter5). 79

PAGE 80

47 ]p.112,[ 8 ]p.267).Let(;B();)beameasurespace.Suppose,t:!CbesuchthatZZjt(x;y)j2d(x)d(y)<1:ThentheassociatedHilbert-SchmidtintegraloperatoristheoperatorT:L2(;C)!L2(;C)givenby(Tf)(x)=Zt(x;y)f(y)d(y):

PAGE 81

PAGE 82

MeynandTweedie 's[ 1993 ]Proposition17.4.1itfollowsthatforHarrisergodicMarkovchainstheonlyfunctionsinL2()satisfyingtheequationlPf=0are-almosteverywhereconstant.SinceL20()doesnotincludefunctionswhichare-almosteverywhereconstant(non-zero)andsinceweconsiderHarrisergodicMarkovchains,theoperatorlPisone-to-oneinourcase.But,lPmightnotbeanontooperatorandsoitmightnotbeinvertible.FromthedenitionofspectrumofanoperatorweknowthatlPisinvertiblei162(P).Since(P)isaclosedsetitfollowsthat162(P)ifandonlyifsup(P)<1.HencetheLaplacianoperator,lP,isinvertibleisup(P)<1.RecallfromChapter2thatwesayaCLTholdsfor 13 ]thataCLTholdsforanyfunctionflyinginR(lP),therangeoflP,andthecorrespondingasymptoticvarianceintheCLTisgivenbyv(f;P)=kgk2kPgk2;wheref=lPgg2L20().So,iflPisinvertible,thenCLTholdsforeveryfunctionf2L20().Also,fromTheorem 7 weknowthatifkPk<1,thenlPisalwaysinvertible.Hence,wehavethefollowingproposition.

PAGE 83

32 ].From[ 13 ]itthenfollowsthatCLTholdsforallf2D(l1P)andbylemma3.2ofMiraandGeyer[ 32 ]wehavev(f;P)=(f;[2l1PI]f);8f2D(l1P):InChapter2,wedenedtheeciencyorderingdueto[ 32 ].TheproblemwitheciencyorderingisthatevenwhenweconcludePEQ,itmighthappenthatv(f;P)=v(f;Q)foreveryf2L2()thatareofinterest.Inthischapter,wediscussconditionsthatguarantee 5{1 )implies( 5{2 ).But,( 5{2 )doesnotnecessarilyimply( 5{1 ).LetlPandlQbetheLaplacianoperatorscorrespondingtotheMarkovoperatorsPandQ.IfweassumethatbothlPandlQareinvertiblethenitfollowsthatD(l1P)=L20()=D(l1Q).Hence,wehavev(f;P)=(f;[2l1PI]f);8f2L20()andv(f;Q)=(f;[2l1QI]f);8f2L20(): 2I0 83

PAGE 84

2I0: 32 ]. 3 ,ifwereplacetheinequalitiesbystrictinequalitiesandcby1,wegetconditionsforstricteciency.FortherestofthissectionweassumethatwehavereversibleMarkovchains.WeknowthattheMarkovtransitionoperator,P,ofareversibleMarkovchainisself-adjoint.So,wecanusespectraltheorytoanalyzereversibleMarkovchains.Let 24 ]showthatifv(f;P)in( 5{3 )isnite,thenthereisaCLTforfwiththeasymptoticvariance,v(f;P)givenby( 5{3 ).Nowonwards,wedenoteEf;f()byEfP().RecallfromChapter2thatPissaidtobebetterthanQintheeciencyorderingwrittenPEQ,ifv(f;P)v(f;Q)foreveryf2L2().ItcanbeeasilyshownthatthedenitionofeciencyorderingisunchangedifthespaceL2()isreplacedwithL20().Itistemptingtoconcludefrom( 5{3 )andthespectraltheoremthatPEQiI+Q IQI+P IP.But,since(P)[1;1],thefunctionh(x)=1+x 84

PAGE 85

32 ]. 35 ],Tierney[ 50 ]).Also,Tierney[ 50 ]showsthatPeskunorderingimpliescovarianceordering.So,P1I,i.e.,lP=IPisapositiveoperator.Fromourdiscussionintheprevioussectionitthenfollowsthatthereexistsauniquesquarerootl1 2PoflP.KipnisandVaradhan[ 24 ]showthatv(f;P)<1ifandonlyiff2R(l1 2P).Notethat,thisdoesnotcontradicttheresultofGordinandLifsic[ 13 ]sinceR(lP)R(l1 2P).MiraandGeyer[ 32 ]provethatP1QifandonlyifPEQ.WegiveanalternativeproofforP1QimpliesPEQ.ThefollowingtheoremisfromBendatandSherman[ 3 ][seealso 32 ,Theorem4.1]. 3 ]mentionthatthefunctionh(x)=ax+b cx+dwithadbc>0satisestheconditionofTheorem 16 eitherinx>d corx
PAGE 86

cx+d,wehaveadbc=2(>0)andd c=1 16 .Since,weareassumingthatP1Qi.e.,QP,byTheorem 16 wehaveh(Q)h(P)for2(0;1):For2(0;1),h(x)isalsoaboundedBorelfunctionon[-1,1].Sowecanusethespectraltheoremtodenetheoperatorh(P).Byspectraltheorem,wehave(h(P)f;f)=Z(P)h()EfP(d):Let+(P)=(P)T[0;1]and(P)=(P)T[1;0)thenwecanwriteZ(P)h()EfP(d)=Z+(P)h()EfP(d)+Z(P)h()EfP(d):Notethat,d dh()=2

PAGE 87

Proof. 16 wehaveQP)I+Q IQI+P IP:Since(P)<1,thefunctionh(x)iscontinuouson(P).ByTheorem 14 ,weknowthath(P)isaboundedoperator.SinceweconsiderthescalareldtobeR,byTheorem 14 ,wealsoknowthath(P)isself-adjoint.Letg(x)=x1 16 ongtoshowI+Q IQI+P IP)QP:Hence,wegetQP,I+Q IQI+P IP:Thefunctionh(x)=1+x IP.Then,by( 5{3 )wehaveI+Q IQI+P IP,PEQ:Hence,P1Q,QP,I+Q IQI+P IP,PEQ:

PAGE 88

2 weknowthatifkPk<1,thenCLTholdsforeverysquareintegrablefunctionf.ItiseasytoseethatifPisacompactoperatorthenkPk<1.So,inordertoestablishCLTforaMarkovchainwecanshowthatthecorrespondingMarkovtransitionoperator,P,iscompact.InSection 5.1 ,wementionedthatanyHilbertSchmidtoperatoriscompact.RecallfromChapter2thatifareversibleMarkovchainisgeometricallyergodic,thentheCLTholdsforeverysquareintegrablefunctionf[ 41 ].AlsowementionedbeforethatareversibleMarkovchain,P,isgeometricallyergodicifandonlyifkPk<1[ 41 44 ].SchervishandCarlin[ 48 ]andLiuetal.[ 26 ])haveprovedgeometricergodicityofcertainMarkovchainsbyestablishingthatthecorrespondingMarkovoperator,P,isHilbertSchmidt. 88

PAGE 89

Herewestate ChenandShao 's[ 2000 ]necessaryandsucientconditionsforc1(y)<1aswellasasimplemethodforcheckingtheseconditions.LetXdenotethenpmatrixwhoseithrowisxTiandletWdenoteannpmatrixwhoseithrowiswTi,wherewi=8><>:xiifyi=0xiifyi=1: 6 ]Thefunctionc1(y)isniteifandonlyif (i) thedesignmatrixXhasfullcolumnrank,and (ii) thereexistsavectora=(a1;:::;an)TwithstrictlypositivecomponentssuchthatWTa=0.AssumingthatXhasfullcolumnrank,thesecondconditionofProposition 5 canbestraightforwardlycheckedwithasimplelinearprogramimplementableintheRprogramminglanguage[ 18 ]usingthe\simplex"functionfromthe\boot"library.Let1andJdenoteacolumnvectorandamatrixof1s,respectively.Thelinearprogramcallsformaximizing1Tasubjectto 5 issatisedandc1(y)<1.Moreover,itisstraightforwardtoshowthatifacontainsoneormorezeros,thentheredoesnotexistanawithallpositiveelementssuchthatWTa=0,soc1(y)=1. 89

PAGE 90

[1] Albert,J.H.andChib,S.(1993),\Bayesiananalysisofbinaryandpolychotomousresponsedata,"JournaloftheAmericanStatisticalAssociation,88,669{679. [2] Arnold,S.F.(1981),TheTheoryofLinearModelsandMultivariateAnalysis,Wiley,NewYork. [3] Bendat,J.andSherman,S.(1955),\Monotoneandconvexoperatorfunctions,"TransactionsoftheAmericanMathematicalSociety,79,58{71. [4] Casella,G.andGeorge,E.(1992),\ExplainingtheGibbssampler,"TheAmericanStatistician,46,167{174. [5] Chan,K.S.andGeyer,C.J.(1994),\Discussionof\MarkovchainsforExploringPosteriorDistributions","TheAnnalsofStatistics,22,1747{1757. [6] Chen,M.-H.andShao,Q.-M.(2000),\Proprietyofposteriordistributionfordichotomousquantalresponsemodels,"ProceedingsoftheAmericanMathemati-calSociety,129,293{302. [7] Chib,S.andGreenberg,E.(1995),\UnderstandingtheMetropolis-Hastingsalgorithm,"TheAmericanStatistician,49,327{335. [8] Conway,J.B.(1990),ACourseinFunctionalAnalysis,Springer-Verlag,NewYork,2nded. [9] Devito,C.L.(1990),FunctionalAnalysisandLinearOperatorTheory,Addison-WesleyPublishingCompany. [10] Eaton,M.L.(1989),GroupInvarianceApplicationsinStatistics,Hayward,CaliforniaandAlexandria,Virginia:InstituteofMathematicalStatisticsandtheAmericanStatisticalAssociation. [11] Feller,W.(1968),AnIntroductiontoProbabilityTheoryanditsApplications,vol.I,NewYork:JohnWiley&Sons,3rded. [12] Fernandez,C.andSteel,M.F.J.(1999),\MultivariateStudent-tregressionmodels:Pitfallsandinference,"Biometrika,86,153{167. [13] Gordin,M.I.andLifsic,B.A.(1978),\ThecentrallimittheoremforstationaryMarkovprocesses,"SovietMathematics.Doklady,19,392{394. [14] Hobert,J.P.(2001),\Discussionof\Theartofdataaugmentation"byD.A.vanDykandX.-L.Meng,"JouranlofComputationalandGraphicalStatistics,10,59{68. [15] Hobert,J.P.andGeyer,C.J.(1998),\GeometricErgodicityofGibbsandBlockGibbsSamplersforaHierarchicalRandomEectsModel,"JournalofMultivariateAnalysis,67,414{430. 90

PAGE 91

Hobert,J.P.,Jones,G.L.,Presnell,B.,andRosenthal,J.S.(2002),\OntheapplicabilityofregenerativesimulationinMarkovchainMonteCarlo,"Biometrika,89,731{743. [17] Hobert,J.P.andMarchev,D.(2008),\Atheoreticalcomparisonofthedataaugmentation,marginalaugmentationandPX-DAalgorithms,"TheAnnalsofStatistics,36,532{554. [18] Ihaka,R.andGentleman,R.(1996),\R:ALanguageforDataAnalysisandGraphics,"JournalofComputationalandGraphicalStatistics,5,299{314. [19] Johnson,N.L.andKotz,S.(1970),ContinuousUnivariateDistributions-1,JohnWiley&Sons. [20] Jones,G.L.(2004),\OntheMarkovchaincentrallimittheorem,"ProbabilitySurveys,1,299{320. [21] Jones,G.L.,Haran,M.,Cao,B.S.,andNeath,R.(2006),\Fixed-widthoutputanalysisforMarkovchainMonteCarlo,"JournaloftheAmericanStatisticalAssocia-tion,101,1537{1547. [22] Jones,G.L.andHobert,J.P.(2001),\HonestexplorationofintractableprobabilitydistributionsviaMarkovchainMonteCarlo,"StatisticalScience,16,312{34. [23] |(2004),\Sucientburn-inforGibbssamplersforahierarchicalrandomeectsmodel,"TheAnnalsofStatistics,32,784{817. [24] Kipnis,C.andVaradhan,S.R.S.(1986),\CentrallimittheoremforadditivefunctionalsofreversibleMarkovprocessesandapplicationstosimpleexclusions,"CommunicationsinMathematicalPhysics,104,1{19. [25] Liu,J.S.,Wong,W.H.,andKong,A.(1994),\CovarianceStructureoftheGibbsSamplerwithApplicationstoComparisonsofEstimatorsandAugmentationSchemes,"Biometrika,81,27{40. [26] |(1995),\CovarianceStructureandConvergenceRateoftheGibbsSamplerwithVariousScans,"JournaloftheRoyalStatisticalSociety,SeriesB,57,157{169. [27] Liu,J.S.andWu,Y.N.(1999),\ParameterExpansionforDataAugmentation,"JournaloftheAmericanStatisticalAssociation,94,1264{1274. [28] Marchev,D.andHobert,J.P.(2004),\GeometricergodicityofvanDykandMeng'salgorithmforthemultivariateStudent'stmodel,"JournaloftheAmericanStatisticalAssociation,99,228{238. [29] Mardia,K.,Kent,J.,andBibby,J.(1979),MultivariateAnalysis,London:Academicpress. 91

PAGE 92

Meyn,S.P.andTweedie,R.L.(1993),MarkovChainsandStochasticStability,London:SpringerVerlag. [31] |(1994),\ComputableboundsforgeometricconvergenceratesofMarkovchains,"TheAnnalsofAppliedProbability,4,981{1011. [32] Mira,A.andGeyer,C.J.(1999),\OrderingMonteCarloMarkovchains,"Tech.Rep.No.632,SchoolofStatistics,UniversityofMinnesota. [33] Mykland,P.,Tierney,L.,andYu,B.(1995),\RegenerationinMarkovchainSamplers,"JournaloftheAmericanStatisticalAssociation,90,233{41. [34] Nummelin,E.(1984),GeneralIrreducibleMarkovChainsandNon-negativeOpera-tors,London:CambridgeUniversityPress. [35] Peskun,P.H.(1973),\OptimumMonteCarlosamplingusingMarkovchains,"Biometrika,60,607{612. [36] RDevelopmentCoreTeam(2006),R:ALanguageandEnvironmentforStatisticalComputing,RFoundationforStatisticalComputing,Vienna,Austria. [37] Retherford,J.R.(1993),HilbertSpace:CompactOperatorsandtheTracetheorem,CambridgeUniversityPress. [38] Robert,C.andCasella,G.(2004),MonteCarloStatisticalMethods,Springer,NewYork. [39] Robert,C.P.(1995),\Simulationoftruncatednormalvariables,"StatisticsandComputing,5,121{125. [40] Roberts,G.andTweedie,R.(1999),\BoundsonregenerationtimesandconvergenceratesforMarkovchains,"StochasticProcessesandtheirApplications,80,221{229.Corrigendum(2001)91:337{338. [41] Roberts,G.O.andRosenthal,J.S.(1997),\GeometricergodicityandhybridMarkovchains,"ElectronicCommunicationsinProbability,2,13{25. [42] |(2001),\MarkovChainsandde-initializingprocesses,"ScandinavianJournalofStatistics,28,489{504. [43] |(2004),\GeneralStateSpaceMarkovChainsandMCMCAlgorithms,"ProbabilitySurveys,1,20{71. [44] Roberts,G.O.andTweedie,R.L.(2001),\GeometricL2andL1convergenceareequivalentforreversibleMarkovchains,"JournalofAppliedProbability,38A,37{41. [45] Rosenthal,J.S.(1995),\MinorizationConditionsandConvergenceRatesforMarkovChainMonteCarlo,"JournaloftheAmericanStatisticalAssociation,90,558{566. 92

PAGE 93

Roy,V.andHobert,J.P.(2007),\ConvergenceratesandasymptoticstandarderrorsforMCMCalgorithmsforBayesianprobitregression,"JournaloftheRoyalStatisticalSociety,SeriesB,69,607{623. [47] Rudin,W.(1991),FunctionalAnalysis,McGraw-Hill,2nded. [48] Schervish,M.J.andCarlin,B.P.(1992),\OntheConvergenceofSuccessiveSubstitutionSampling,"JouranlofComputationalandGraphicalStatistics,1,111{127. [49] Tierney,L.(1994),\Markovchainsforexploringposteriordistributions(withdiscussion),"TheAnnalsofStatistics,22,1701{1762. [50] |(1998),\ANoteonMetropolis-HastingsKernelsforGeneralStateSpaces,"TheAnnalsofAppliedProbability,8,1{9. [51] vanDyk,D.A.andMeng,X.-L.(2001),\TheArtofDataAugmentation(withDiscussion),"JournalofComputationalandGraphicalStatistics,10,1{50. 93

PAGE 94