メトロポリス・ヘイスティング・アルゴリズム

Bayesian Biostatistics by E Lesaffre & AB Lawson, Chapter 6 学習ノート
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
とりあえず、プログラム実行してから考えることとする。

> # MH SAMPLING with NONINFORMATIVE PRIOR
> 
> 
> lposterior <- function(nsample, samplemean, samplevar, mu,sigma2)
+ {
+     term1 <- -((nsample-1)*samplevar+nsample*(samplemean-mu)^2)/(2*sigma2)
+     term2 <- -0.5*(nsample+2)*log(sigma2)
+     lposterior <- term1+term2
+ }
> 
> attach(ALP)
> 
> objects(2)
[1] "alkfos"  "artikel"
> alp <- alkfos[artikel==0]
> talp <- 100*alp^(-1/2)
> 
> mtalp <- mean(talp)
> vartalp <- var(talp)
> sdtalp <- sqrt(var(talp))
> ntalp <- length(talp)
> 
> ntalp;mtalp;sdtalp;vartalp
[1] 250
[1] 7.108724
[1] 1.365344
[1] 1.864165

> # MH SAMPLING with NONINFORMATIVE PRIOR

> lposterior <- function(nsample, samplemean, samplevar, mu,sigma2)

+ {

+ term1 <- -((nsample-1)*samplevar+nsample*(samplemean-mu)^2)/(2*sigma2)

+ term2 <- -0.5*(nsample+2)*log(sigma2)

+ lposterior <- term1+term2

+ }

> attach(ALP)

> objects(2)

[1] "alkfos" "artikel"

> alp <- alkfos[artikel==0]

> talp <- 100*alp^(-1/2)

> mtalp <- mean(talp)

> vartalp <- var(talp)

> sdtalp <- sqrt(var(talp))

> ntalp <- length(talp)

> ntalp;mtalp;sdtalp;vartalp

[1] 250

[1] 7.108724

[1] 1.365344

[1] 1.864165

> # MH sampling of normal data set
> 
> nchain <- 1500
> 
> theta.1 <- rep(0,nchain)
> theta.2 <- rep(0,nchain)
> 
> theta.accept <- rep(0,nchain)
> 
> # Figure 6.11
> # Use the sigma and tau values equal to 0.03
> 
> sigma2s <- 0.03
> tau2s <- 0.03
> set.seed(729)
> 
> # Figure 6.14
> # Use the sigma and tau values equal to 0.001
> # Thus put above three lines in comment and rerun all programs
> 
> sigma2s <- 0.001
> tau2s <- 0.001
> set.seed(729)
> 
> ichain <- 1
> accept <- 0
> theta.1[ichain] <- 6.5
> theta.2[ichain] <- 2
> 
> theta.accept[ichain] <- 1
> 
> 
> par(mfrow=c(1,1))
> par(pty="s")
> par(mar=c(5,5,2,4))
> 
> plot(theta.1[ichain],theta.2[ichain],xlab=expression(paste(mu)),
+      ylab=expression(sigma^{2}),type="n",bty="o",
+      xlim=c(6.5,7.5),ylim=c(1.4,2.65),main="",cex.lab=1.5,cex.axis=1.3)
> 
> points(theta.1[ichain],theta.2[ichain],cex=1.2,pch=15,col="blue")
> text(6.55,2.6,"(a)",cex=1.3)
>

> # MH sampling of normal data set

> nchain <- 1500

> theta.1 <- rep(0,nchain)

> theta.2 <- rep(0,nchain)

> theta.accept <- rep(0,nchain)

> # Figure 6.11

> # Use the sigma and tau values equal to 0.03

> sigma2s <- 0.03

> tau2s <- 0.03

> set.seed(729)

> # Figure 6.14

> # Use the sigma and tau values equal to 0.001

> # Thus put above three lines in comment and rerun all programs

> sigma2s <- 0.001

> tau2s <- 0.001

> set.seed(729)

> ichain <- 1

> accept <- 0

> theta.1[ichain] <- 6.5

> theta.2[ichain] <- 2

> theta.accept[ichain] <- 1

> par(mfrow=c(1,1))

> par(pty="s")

> par(mar=c(5,5,2,4))

> plot(theta.1[ichain],theta.2[ichain],xlab=expression(paste(mu)),

+ ylab=expression(sigma^{2}),type="n",bty="o",

+ xlim=c(6.5,7.5),ylim=c(1.4,2.65),main="",cex.lab=1.5,cex.axis=1.3)

> points(theta.1[ichain],theta.2[ichain],cex=1.2,pch=15,col="blue")

> text(6.55,2.6,"(a)",cex=1.3)

> while (ichain < nchain) {   
+     theta.1.tilde <- rnorm(1,theta.1[ichain],sqrt(sigma2s))
+     theta.2.tilde <- rnorm(1,theta.2[ichain],sqrt(tau2s))
+     
+     lpostnew <- lposterior(ntalp, mtalp, vartalp, theta.1.tilde,   theta.2.tilde)
+     lpostold <- lposterior(ntalp, mtalp, vartalp, theta.1[ichain], theta.2[ichain])
+     r <- exp(lpostnew-lpostold)
+     
+     u <- runif(1,0,1)
+     
+     if(u < r) {theta.1[ichain+1] <- theta.1.tilde; 
+     theta.2[ichain+1] <- theta.2.tilde; 
+     accept=accept+1; 
+     theta.accept[ichain+1] <- 1}
+     else {theta.1[ichain+1] <- theta.1[ichain]; 
+     theta.2[ichain+1] <- theta.2[ichain]; 
+     theta.accept[ichain+1] <- 0}
+     
+     if (ichain <= 15){
+         segments(theta.1[ichain],theta.2[ichain],theta.1[ichain+1],theta.2[ichain+1])
+         
+         # red points are suggested moves, but were not accepted
+         # these points are not displayed in Figures 6.11 and 6.14
+         
+         points(theta.1.tilde,theta.2.tilde,cex=1.8,pch=20,col="red") 
+         points(theta.1[ichain+1],theta.2[ichain+1],cex=1.8,pch=15,col="blue")
+         print(c("theta.1.tilde=", theta.1.tilde))
+         print(c("theta.2.tilde=", theta.2.tilde))
+         print(c("theta1[ichain]=", theta.1[ichain]))
+         print(c("theta2[ichain]=", theta.2[ichain]))
+         print(c("r=",r))
+     }
+     
+     ichain <- ichain +1   
+     
+ }
[1] "theta.1.tilde="   "6.47555478090042"
[1] "theta.2.tilde="   "1.98995455325182"
[1] "theta1[ichain]=" "6.5"            
[1] "theta2[ichain]=" "2"              
[1] "r="                "0.138730006426686"
[1] "theta.1.tilde="   "6.55359799612869"
[1] "theta.2.tilde="  "2.0756284755932"
[1] "theta1[ichain]=" "6.5"            
[1] "theta2[ichain]=" "2"              
[1] "r="               "63.5638743254562"
[1] "theta.1.tilde="   "6.51826420865378"
[1] "theta.2.tilde="   "2.05622829464034"
[1] "theta1[ichain]="  "6.55359799612869"
[1] "theta2[ichain]=" "2.0756284755932"
[1] "r="                 "0.0814690245530072"
[1] "theta.1.tilde="   "6.55679897666421"
[1] "theta.2.tilde="   "2.06009971600029"
[1] "theta1[ichain]="  "6.51826420865378"
[1] "theta2[ichain]="  "2.05622829464034"
[1] "r="               "14.6730553950733"
[1] "theta.1.tilde="   "6.56555412837179"
[1] "theta.2.tilde="   "2.05161733250196"
[1] "theta1[ichain]="  "6.55679897666421"
[1] "theta2[ichain]="  "2.06009971600029"
[1] "r="               "1.75386431547291"
[1] "theta.1.tilde="   "6.53331664995836"
[1] "theta.2.tilde="   "2.00445949942664"
[1] "theta1[ichain]="  "6.56555412837179"
[1] "theta2[ichain]="  "2.05161733250196"
[1] "r="                 "0.0904402175581735"
[1] "theta.1.tilde="   "6.53723882032352"
[1] "theta.2.tilde="   "2.03451061074855"
[1] "theta1[ichain]="  "6.56555412837179"
[1] "theta2[ichain]="  "2.05161733250196"
[1] "r="                "0.137179511182405"
[1] "theta.1.tilde="   "6.57539178365308"
[1] "theta.2.tilde="   "2.08247701839794"
[1] "theta1[ichain]="  "6.53723882032352"
[1] "theta2[ichain]="  "2.03451061074855"
[1] "r="               "14.6398686576767"
[1] "theta.1.tilde="   "6.57486014348542"
[1] "theta.2.tilde="   "2.03475471794679"
[1] "theta1[ichain]="  "6.57539178365308"
[1] "theta2[ichain]="  "2.08247701839794"
[1] "r="                "0.879744590424981"
[1] "theta.1.tilde="   "6.62128168415793"
[1] "theta.2.tilde="   "2.01276516714005"
[1] "theta1[ichain]="  "6.57486014348542"
[1] "theta2[ichain]="  "2.03475471794679"
[1] "r="               "17.7438633629861"
[1] "theta.1.tilde="   "6.63516899348182"
[1] "theta.2.tilde="   "1.99864981040392"
[1] "theta1[ichain]="  "6.62128168415793"
[1] "theta2[ichain]="  "2.01276516714005"
[1] "r="               "2.23186427515578"
[1] "theta.1.tilde="   "6.63364844624157"
[1] "theta.2.tilde="   "2.05080212497638"
[1] "theta1[ichain]="  "6.63516899348182"
[1] "theta2[ichain]="  "1.99864981040392"
[1] "r="                "0.976402020691009"
[1] "theta.1.tilde="   "6.63313903461943"
[1] "theta.2.tilde="   "1.97996807104503"
[1] "theta1[ichain]="  "6.63364844624157"
[1] "theta2[ichain]="  "2.05080212497638"
[1] "r="                "0.867214862858916"
[1] "theta.1.tilde="   "6.65059159697973"
[1] "theta.2.tilde="   "1.92250721992873"
[1] "theta1[ichain]="  "6.63313903461943"
[1] "theta2[ichain]="  "1.97996807104503"
[1] "r="              "2.3164588832004"
[1] "theta.1.tilde="  "6.7156321531923"
[1] "theta.2.tilde="   "1.90464656615876"
[1] "theta1[ichain]="  "6.65059159697973"
[1] "theta2[ichain]="  "1.92250721992873"
[1] "r="               "34.7961839504696"

100

101

102

103

104

105

106

107

108

109

110

111

> while (ichain < nchain) {

+ theta.1.tilde <- rnorm(1,theta.1[ichain],sqrt(sigma2s))

+ theta.2.tilde <- rnorm(1,theta.2[ichain],sqrt(tau2s))

+ lpostnew <- lposterior(ntalp, mtalp, vartalp, theta.1.tilde, theta.2.tilde)

+ lpostold <- lposterior(ntalp, mtalp, vartalp, theta.1[ichain], theta.2[ichain])

+ r <- exp(lpostnew-lpostold)

+ u <- runif(1,0,1)

+ if(u < r) {theta.1[ichain+1] <- theta.1.tilde;

+ theta.2[ichain+1] <- theta.2.tilde;

+ accept=accept+1;

+ theta.accept[ichain+1] <- 1}

+ else {theta.1[ichain+1] <- theta.1[ichain];

+ theta.2[ichain+1] <- theta.2[ichain];

+ theta.accept[ichain+1] <- 0}

+ if (ichain <= 15){

+ segments(theta.1[ichain],theta.2[ichain],theta.1[ichain+1],theta.2[ichain+1])

+ # red points are suggested moves, but were not accepted

+ # these points are not displayed in Figures 6.11 and 6.14

+ points(theta.1.tilde,theta.2.tilde,cex=1.8,pch=20,col="red")

+ points(theta.1[ichain+1],theta.2[ichain+1],cex=1.8,pch=15,col="blue")

+ print(c("theta.1.tilde=", theta.1.tilde))

+ print(c("theta.2.tilde=", theta.2.tilde))

+ print(c("theta1[ichain]=", theta.1[ichain]))

+ print(c("theta2[ichain]=", theta.2[ichain]))

+ print(c("r=",r))

+ }

+ ichain <- ichain +1

+ }

[1] "theta.1.tilde=" "6.47555478090042"

[1] "theta.2.tilde=" "1.98995455325182"

[1] "theta1[ichain]=" "6.5"

[1] "theta2[ichain]=" "2"

[1] "r=" "0.138730006426686"

[1] "theta.1.tilde=" "6.55359799612869"

[1] "theta.2.tilde=" "2.0756284755932"

[1] "theta1[ichain]=" "6.5"

[1] "theta2[ichain]=" "2"

[1] "r=" "63.5638743254562"

[1] "theta.1.tilde=" "6.51826420865378"

[1] "theta.2.tilde=" "2.05622829464034"

[1] "theta1[ichain]=" "6.55359799612869"

[1] "theta2[ichain]=" "2.0756284755932"

[1] "r=" "0.0814690245530072"

[1] "theta.1.tilde=" "6.55679897666421"

[1] "theta.2.tilde=" "2.06009971600029"

[1] "theta1[ichain]=" "6.51826420865378"

[1] "theta2[ichain]=" "2.05622829464034"

[1] "r=" "14.6730553950733"

[1] "theta.1.tilde=" "6.56555412837179"

[1] "theta.2.tilde=" "2.05161733250196"

[1] "theta1[ichain]=" "6.55679897666421"

[1] "theta2[ichain]=" "2.06009971600029"

[1] "r=" "1.75386431547291"

[1] "theta.1.tilde=" "6.53331664995836"

[1] "theta.2.tilde=" "2.00445949942664"

[1] "theta1[ichain]=" "6.56555412837179"

[1] "theta2[ichain]=" "2.05161733250196"

[1] "r=" "0.0904402175581735"

[1] "theta.1.tilde=" "6.53723882032352"

[1] "theta.2.tilde=" "2.03451061074855"

[1] "theta1[ichain]=" "6.56555412837179"

[1] "theta2[ichain]=" "2.05161733250196"

[1] "r=" "0.137179511182405"

[1] "theta.1.tilde=" "6.57539178365308"

[1] "theta.2.tilde=" "2.08247701839794"

[1] "theta1[ichain]=" "6.53723882032352"

[1] "theta2[ichain]=" "2.03451061074855"

[1] "r=" "14.6398686576767"

[1] "theta.1.tilde=" "6.57486014348542"

[1] "theta.2.tilde=" "2.03475471794679"

[1] "theta1[ichain]=" "6.57539178365308"

[1] "theta2[ichain]=" "2.08247701839794"

[1] "r=" "0.879744590424981"

[1] "theta.1.tilde=" "6.62128168415793"

[1] "theta.2.tilde=" "2.01276516714005"

[1] "theta1[ichain]=" "6.57486014348542"

[1] "theta2[ichain]=" "2.03475471794679"

[1] "r=" "17.7438633629861"

[1] "theta.1.tilde=" "6.63516899348182"

[1] "theta.2.tilde=" "1.99864981040392"

[1] "theta1[ichain]=" "6.62128168415793"

[1] "theta2[ichain]=" "2.01276516714005"

[1] "r=" "2.23186427515578"

[1] "theta.1.tilde=" "6.63364844624157"

[1] "theta.2.tilde=" "2.05080212497638"

[1] "theta1[ichain]=" "6.63516899348182"

[1] "theta2[ichain]=" "1.99864981040392"

[1] "r=" "0.976402020691009"

[1] "theta.1.tilde=" "6.63313903461943"

[1] "theta.2.tilde=" "1.97996807104503"

[1] "theta1[ichain]=" "6.63364844624157"

[1] "theta2[ichain]=" "2.05080212497638"

[1] "r=" "0.867214862858916"

[1] "theta.1.tilde=" "6.65059159697973"

[1] "theta.2.tilde=" "1.92250721992873"

[1] "theta1[ichain]=" "6.63313903461943"

[1] "theta2[ichain]=" "1.97996807104503"

[1] "r=" "2.3164588832004"

[1] "theta.1.tilde=" "6.7156321531923"

[1] "theta.2.tilde=" "1.90464656615876"

[1] "theta1[ichain]=" "6.65059159697973"

[1] "theta2[ichain]=" "1.92250721992873"

[1] "r=" "34.7961839504696"

> # Compute acceptance rate
> 
> 
> rate.accept <- accept/nchain
> rate.accept
[1] 0.8326667
> 
> # Show scatterplot of last 1000 sampled values
> 
> nstart <- 501
> 
> burnin <- seq(1,nstart-1)
> 
> 
> 
> plot(theta.1[nstart:nchain],theta.2[nstart:nchain],
+      xlab=expression(mu),
+      ylab=expression(sigma^{2}),type="n",bty="o",
+      xlim=c(6.5,7.5),ylim=c(1.4,2.65),main="",cex.lab=1.5,cex.axis=1.3)
> points(theta.1[nstart:nchain],theta.2[nstart:nchain],cex=0.7,pch=1,col="blue")
> text(6.55,2.6,"(b)",cex=1.3)
> 
>

> # Compute acceptance rate

> rate.accept <- accept/nchain

> rate.accept

[1] 0.8326667

> # Show scatterplot of last 1000 sampled values

> nstart <- 501

> burnin <- seq(1,nstart-1)

> plot(theta.1[nstart:nchain],theta.2[nstart:nchain],

+ xlab=expression(mu),

+ ylab=expression(sigma^{2}),type="n",bty="o",

+ xlim=c(6.5,7.5),ylim=c(1.4,2.65),main="",cex.lab=1.5,cex.axis=1.3)

> points(theta.1[nstart:nchain],theta.2[nstart:nchain],cex=0.7,pch=1,col="blue")

> text(6.55,2.6,"(b)",cex=1.3)

> # Show unique values in the chain
> # Excluding when the chain stayed at the same position
> # This is used to indicate in the trace plots below when the chain did not move
> # Rejections are indicated by red horizontal lines
> 
> unique(theta.1[nstart:nchain])
  [1] 7.105710 7.097377 7.146154 7.126149 7.138248 7.136730 7.193407 7.196787
  [9] 7.161477 7.167565 7.160029 7.149714 7.166578 7.148522 7.148769 7.111158
 [17] 7.066738 7.074376 7.079513 7.105338 7.110835 7.108151 7.077600 7.032159
 [25] 7.030563 7.039744 7.001758 7.011440 7.038886 7.043371 6.984262 7.000522
 [33] 7.000138 7.052493 6.998251 6.931824 6.975568 6.973657 6.986464 7.008727
 [41] 7.023473 6.981957 6.997285 7.026417 7.028393 7.090305 7.195698 7.220242
 [49] 7.193362 7.166483 7.133860 7.160466 7.201392 7.184865 7.175697 7.113790
 [57] 7.095223 7.141158 7.150247 7.168451 7.163010 7.185438 7.125459 7.118104
 [65] 7.070248 7.077363 7.112081 7.137304 7.150435 7.172319 7.162484 7.143403
 [73] 7.150210 7.122165 7.123491 7.091077 7.126873 7.105779 7.139805 7.094731
 [81] 7.096497 7.102668 7.110727 7.167064 7.177816 7.147048 7.137519 7.157454
 [89] 7.130782 7.161615 7.158419 7.078840 7.090250 7.074787 7.050780 7.055069
 [97] 7.103408 7.094454 7.069294 7.070464 7.079295 7.039022 6.995843 6.988863
[105] 6.982750 7.009520 7.019409 7.024520 7.026422 7.055173 7.046609 7.021164
[113] 7.011415 7.086306 7.071017 7.096240 7.118657 7.152039 7.168983 7.161043
[121] 7.134903 7.146904 7.129693 7.145513 7.164515 7.146554 7.179837 7.171510
[129] 7.165563 7.192714 7.195835 7.173599 7.173662 7.210244 7.221889 7.225565
[137] 7.259889 7.278050 7.299470 7.310388 7.351809 7.343904 7.309893 7.321255
[145] 7.336674 7.284471 7.280756 7.255376 7.268188 7.213821 7.199986 7.184075
[153] 7.190026 7.179515 7.206078 7.203055 7.190790 7.204960 7.169081 7.197101
[161] 7.174371 7.167990 7.174018 7.210532 7.249923 7.236865 7.236578 7.196585
[169] 7.185049 7.157782 7.152770 7.136116 7.071804 7.053627 7.021525 7.062086
[177] 7.050403 7.018537 7.038144 7.088847 7.116697 7.148690 7.148674 7.132627
[185] 7.123776 7.131811 7.169534 7.114034 7.103459 7.137737 7.113494 7.151167
[193] 7.118326 7.162041 7.124493 7.093596 7.103819 7.059961 7.126471 7.169115
[201] 7.176733 7.159343 7.172765 7.134412 7.078622 7.092632 7.104736 7.098675
[209] 7.111029 7.118260 7.154540 7.141623 7.137082 7.114960 7.136295 7.190459
[217] 7.203320 7.214676 7.163769 7.156705 7.146943 7.112468 7.107483 7.124075
[225] 7.081330 7.090805 7.102128 7.153387 7.156471 7.146052 7.091781 7.060955
[233] 7.074721 7.101063 7.127645 7.093060 7.124284 7.134584 7.113500 7.127588
[241] 7.081936 7.084418 7.073552 7.048281 7.037618 7.016605 7.081087 7.099355
[249] 7.109486 7.051725 7.069731 7.066463 7.063431 7.082629 7.072206 7.060580
[257] 7.068385 7.062276 7.071833 7.108285 7.103986 7.060831 7.000828 7.007190
[265] 7.010223 7.041649 7.083511 7.075521 7.087726 7.037997 7.012511 7.002025
[273] 7.018571 7.017577 7.075386 7.107743 7.132775 7.108093 7.121088 7.097991
[281] 7.085211 7.079735 7.121678 7.085073 7.095945 7.097626 7.052727 7.030969
[289] 7.023504 7.065117 7.066804 7.121666 7.170794 7.135517 7.131297 7.125362
[297] 7.101836 7.111131 7.066778 7.118803 7.108834 7.118440 7.085606 7.140762
[305] 7.157767 7.198726 7.217093 7.144058 7.089210 7.103224 7.116330 7.133247
[313] 7.182631 7.219941 7.234179 7.233300 7.227361 7.293898 7.297476 7.300428
[321] 7.289837 7.273858 7.288363 7.244932 7.237168 7.271674 7.255755 7.257549
[329] 7.274788 7.249515 7.244348 7.247958 7.298818 7.295422 7.300232 7.347096
[337] 7.305467 7.293871 7.279378 7.288125 7.257229 7.262147 7.230256 7.170210
[345] 7.137925 7.104712 7.121480 7.118843 7.117378 7.090968 7.084255 7.108023
[353] 7.051757 7.028121 7.054859 7.078414 7.114706 7.079548 7.114003 7.078881
[361] 7.090100 7.129850 7.126992 7.139105 7.094334 7.087728 7.105869 7.068680
[369] 7.077331 7.074386 7.037065 7.015845 7.031067 7.045745 7.045352 7.067747
[377] 7.079858 7.077529 7.085807 7.101523 7.146237 7.158741 7.093481 7.064842
[385] 7.096556 7.130847 7.141951 7.115001 7.112850 7.150857 7.145183 7.094028
[393] 7.118288 7.150233 7.145857 7.144554 7.146988 7.146072 7.176767 7.107342
[401] 7.123751 7.149823 7.156241 7.111539 7.146516 7.174778 7.180236 7.151633
[409] 7.177896 7.097421 7.134737 7.139064 7.180692 7.176342 7.181836 7.204021
[417] 7.208299 7.146985 7.148830 7.107342 7.088155 7.052253 7.006641 7.015996
[425] 7.046809 7.061120 7.032716 7.041052 7.064903 7.059522 7.091234 7.090501
[433] 7.081098 7.045878 7.059936 7.085968 7.130180 7.102629 7.103868 7.103656
[441] 7.117224 7.131319 7.098237 7.118564 7.101916 7.079309 7.078102 7.027398
[449] 7.017138 7.023637 7.091488 7.122048 7.133928 7.131692 7.159420 7.169551
[457] 7.146890 7.142987 7.140642 7.126635 7.175829 7.154942 7.128365 7.181859
[465] 7.175072 7.169410 7.238363 7.237141 7.289526 7.306859 7.273925 7.263208
[473] 7.243585 7.268352 7.224143 7.250465 7.270233 7.255779 7.240344 7.291077
[481] 7.307398 7.266715 7.266844 7.284025 7.303616 7.278101 7.329640 7.312715
[489] 7.345105 7.305630 7.276796 7.269477 7.256781 7.295397 7.296595 7.320085
[497] 7.331981 7.331155 7.333449 7.307147 7.316371 7.328005 7.300980 7.281100
[505] 7.288726 7.225282 7.232577 7.229973 7.249190 7.323170 7.325466 7.359743
[513] 7.354087 7.314519 7.316979 7.316903 7.295794 7.266593 7.290003 7.309931
[521] 7.269465 7.240654 7.244580 7.236161 7.220369 7.216624 7.214971 7.182581
[529] 7.164027 7.127042 7.135585 7.153463 7.168115 7.128040 7.155798 7.135208
[537] 7.177279 7.221284 7.256975 7.279565 7.272557 7.250464 7.257887 7.278663
[545] 7.233752 7.227740 7.232923 7.291705 7.251748 7.216629 7.246354 7.246186
[553] 7.222686 7.240670 7.240207 7.191035 7.232885 7.270038 7.227842 7.217455
[561] 7.170541 7.178615 7.206536 7.192185 7.185152 7.186717 7.157135 7.174848
[569] 7.183865 7.156051 7.142732 7.105314 7.068098 7.087493 7.020719 7.031725
[577] 7.016702 7.015673 6.999963 7.045578 7.069696 7.110356 7.116568 7.111650
[585] 7.046876 7.098489 7.098522 7.104072 7.135575 7.074844 7.090479 7.100506
[593] 7.121507 7.123277 7.090857 7.085827 7.079468 7.036278 7.043417 7.009129
[601] 7.006092 6.995034 6.997751 6.961613 6.962277 7.000038 6.992817 7.028478
[609] 7.030588 7.004229 6.975593 6.957772 6.948820 6.938782 6.912071 6.900986
[617] 6.913416 6.942149 6.958355 7.004749 6.981839 7.010009 7.034225 7.090857
[625] 7.087262 7.141156 7.165744 7.170199 7.123199 7.155267 7.092919 7.153585
[633] 7.159173 7.095379 7.123217 7.091858 7.103740 7.071565 7.035169 7.066140
[641] 7.095298 7.084756 7.093266 7.078597 7.089700 7.139309 7.117069 7.137335
[649] 7.136439 7.116762 7.173263 7.133395 7.147188 7.203146 7.242657 7.250927
[657] 7.201513 7.188890 7.169367 7.157233 7.103280 7.040244 7.070839 7.083728
[665] 7.058956 7.024028 7.031458 7.069782 7.106274 7.130350 7.126406 7.142285
[673] 7.132842 7.136819 7.126732 7.165840 7.160837 7.147163 7.121239 7.107829
[681] 7.094635 7.127397 7.129823 7.144622 7.129571 7.127359 7.095176 7.072701
[689] 7.025348 7.026530 7.053881 7.038662 7.077207 7.065855 7.092630 7.090197
[697] 7.104571 7.065959 7.086947 7.086198 7.061146 7.071278 7.087874 7.075494
[705] 7.084810 7.071096 7.017197 7.028904 7.045269 6.995492 6.993004 7.002403
[713] 7.070153 7.048646 7.053306 7.068167 7.065019 7.051720 7.029827 7.067742
[721] 7.101737 7.074833 7.109695 7.117952 7.171604 7.157373 7.122277 7.102073
[729] 7.122280 7.092120 7.065516 7.036670 7.006340 7.011977 7.011673 6.996891
[737] 7.014798 7.015064 6.995969 6.958202 6.973073 6.979068 6.937016 6.969629
[745] 6.984088 6.976605 6.951565 6.942892 6.924040 6.907088 6.912183 6.905260
[753] 6.903310 6.956582 6.930251 6.886287 6.911363 6.924639 6.914682 6.911524
[761] 6.893737 6.899142 6.907357 6.914254 6.927394 6.909583 6.892624 6.903670
[769] 6.903444 6.902892 6.918966 6.955622 6.950996 6.933009 6.895144 6.931373
[777] 6.889020 6.882409 6.874159 6.888222 6.897185 6.883539 6.915821 6.932445
[785] 6.940479 6.929080 6.931299 6.964280 6.999782 6.989388 6.981434 6.968265
[793] 7.003788 6.995484 7.003156 6.954861 6.982269 6.972628 6.998854 7.000844
[801] 6.959610 6.940751 6.952307 6.919167 6.934823 6.991495 7.009198 7.024987
[809] 7.010942 7.009779 7.026069 7.053679 7.069666 7.098924 7.100481 7.088994
[817] 7.088969 7.037677 7.023283 7.052647 7.074681 7.051656 7.132211 7.057101
[825] 7.065356 7.075734 7.102188 7.155825 7.167115 7.153556 7.088877 7.113851
[833] 7.099440 7.106808 7.036434 7.001959 7.003078 7.086281 7.089996 7.053244
[841] 7.093258
> 
> theta.accept
   [1] 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 0 0 0 1 1 1 0 1 1 1 1
  [37] 0 1 0 1 1 1 0 0 1 1 1 1 0 0 1 0 1 1 0 1 0 1 0 1 1 0 1 0 1 1 1 1 0 1 1 1
  [73] 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
 [109] 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1
 [145] 0 1 1 0 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
 [181] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1
 [217] 1 1 1 1 1 0 1 1 0 0 0 0 1 1 1 0 1 1 0 0 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1
 [253] 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [289] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 1 1
 [325] 1 0 1 1 1 1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1
 [361] 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1
 [397] 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 0
 [433] 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 0
 [469] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1
 [505] 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1
 [541] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1
 [577] 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 1 1 0
 [613] 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [649] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [685] 1 1 1 0 0 1 0 1 1 1 0 1 1 1 0 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
 [721] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1
 [757] 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 1 1 1 1 1
 [793] 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 0 0 1 0 0 1 1 0 1 1 1 0 1 1
 [829] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1
 [865] 0 1 1 1 1 1 1 0 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1
 [901] 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 1 1 0
 [937] 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 1 1 1
 [973] 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
 [ reached getOption("max.print") -- omitted 500 entries ]
> 
> colchain <- rep("blue",nchain)
> 
> for (ichain in 1: nchain){
+     if (theta.accept[ichain] == 0) colchain[ichain]="red"
+ }
> 
> # Traceplots
> # Figure 6.13
> 
> par(mar=c(4,5,4,2))
> 
> par(mfrow=c(1,1))
> par(pty="m")
> 
> 
> plot(nstart:nchain,theta.1[nstart:nchain],xlab="Iteration",
+      ylab=expression(paste(mu)),cex.lab=1.5,cex=1.3,cex.axis=1.3,
+      type="n",main="",bty="o")
> for (ichain in nstart:(nchain-1)){
+     segments(ichain,theta.1[ichain],ichain+1,theta.1[ichain+1],
+              col=colchain[ichain+1],lwd=2)
+ }
> text(500,7.3,"(a)",cex=1.5)
> 
> 
> plot(nstart:nchain,theta.2[nstart:nchain],xlab="Iteration",
+      ylab=expression(paste(sigma^{2})),cex.lab=1.5,cex=1.5,cex.axis=1.3,
+      type="n",main="",bty="o")
> for (ichain in nstart:(nchain-1)){
+     segments(ichain,theta.2[ichain],ichain+1,theta.2[ichain+1],
+              col=colchain[ichain+1],lwd=2)
+ }
> text(500,2.37,"(b)",cex=1.5)
> 
>

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

> # Show unique values in the chain

> # Excluding when the chain stayed at the same position

> # This is used to indicate in the trace plots below when the chain did not move

> # Rejections are indicated by red horizontal lines

> unique(theta.1[nstart:nchain])

[1] 7.105710 7.097377 7.146154 7.126149 7.138248 7.136730 7.193407 7.196787

[9] 7.161477 7.167565 7.160029 7.149714 7.166578 7.148522 7.148769 7.111158

[17] 7.066738 7.074376 7.079513 7.105338 7.110835 7.108151 7.077600 7.032159

[25] 7.030563 7.039744 7.001758 7.011440 7.038886 7.043371 6.984262 7.000522

[33] 7.000138 7.052493 6.998251 6.931824 6.975568 6.973657 6.986464 7.008727

[41] 7.023473 6.981957 6.997285 7.026417 7.028393 7.090305 7.195698 7.220242

[49] 7.193362 7.166483 7.133860 7.160466 7.201392 7.184865 7.175697 7.113790

[57] 7.095223 7.141158 7.150247 7.168451 7.163010 7.185438 7.125459 7.118104

[65] 7.070248 7.077363 7.112081 7.137304 7.150435 7.172319 7.162484 7.143403

[73] 7.150210 7.122165 7.123491 7.091077 7.126873 7.105779 7.139805 7.094731

[81] 7.096497 7.102668 7.110727 7.167064 7.177816 7.147048 7.137519 7.157454

[89] 7.130782 7.161615 7.158419 7.078840 7.090250 7.074787 7.050780 7.055069

[97] 7.103408 7.094454 7.069294 7.070464 7.079295 7.039022 6.995843 6.988863

[105] 6.982750 7.009520 7.019409 7.024520 7.026422 7.055173 7.046609 7.021164

[113] 7.011415 7.086306 7.071017 7.096240 7.118657 7.152039 7.168983 7.161043

[121] 7.134903 7.146904 7.129693 7.145513 7.164515 7.146554 7.179837 7.171510

[129] 7.165563 7.192714 7.195835 7.173599 7.173662 7.210244 7.221889 7.225565

[137] 7.259889 7.278050 7.299470 7.310388 7.351809 7.343904 7.309893 7.321255

[145] 7.336674 7.284471 7.280756 7.255376 7.268188 7.213821 7.199986 7.184075

[153] 7.190026 7.179515 7.206078 7.203055 7.190790 7.204960 7.169081 7.197101

[161] 7.174371 7.167990 7.174018 7.210532 7.249923 7.236865 7.236578 7.196585

[169] 7.185049 7.157782 7.152770 7.136116 7.071804 7.053627 7.021525 7.062086

[177] 7.050403 7.018537 7.038144 7.088847 7.116697 7.148690 7.148674 7.132627

[185] 7.123776 7.131811 7.169534 7.114034 7.103459 7.137737 7.113494 7.151167

[193] 7.118326 7.162041 7.124493 7.093596 7.103819 7.059961 7.126471 7.169115

[201] 7.176733 7.159343 7.172765 7.134412 7.078622 7.092632 7.104736 7.098675

[209] 7.111029 7.118260 7.154540 7.141623 7.137082 7.114960 7.136295 7.190459

[217] 7.203320 7.214676 7.163769 7.156705 7.146943 7.112468 7.107483 7.124075

[225] 7.081330 7.090805 7.102128 7.153387 7.156471 7.146052 7.091781 7.060955

[233] 7.074721 7.101063 7.127645 7.093060 7.124284 7.134584 7.113500 7.127588

[241] 7.081936 7.084418 7.073552 7.048281 7.037618 7.016605 7.081087 7.099355

[249] 7.109486 7.051725 7.069731 7.066463 7.063431 7.082629 7.072206 7.060580

[257] 7.068385 7.062276 7.071833 7.108285 7.103986 7.060831 7.000828 7.007190

[265] 7.010223 7.041649 7.083511 7.075521 7.087726 7.037997 7.012511 7.002025

[273] 7.018571 7.017577 7.075386 7.107743 7.132775 7.108093 7.121088 7.097991

[281] 7.085211 7.079735 7.121678 7.085073 7.095945 7.097626 7.052727 7.030969

[289] 7.023504 7.065117 7.066804 7.121666 7.170794 7.135517 7.131297 7.125362

[297] 7.101836 7.111131 7.066778 7.118803 7.108834 7.118440 7.085606 7.140762

[305] 7.157767 7.198726 7.217093 7.144058 7.089210 7.103224 7.116330 7.133247

[313] 7.182631 7.219941 7.234179 7.233300 7.227361 7.293898 7.297476 7.300428

[321] 7.289837 7.273858 7.288363 7.244932 7.237168 7.271674 7.255755 7.257549

[329] 7.274788 7.249515 7.244348 7.247958 7.298818 7.295422 7.300232 7.347096

[337] 7.305467 7.293871 7.279378 7.288125 7.257229 7.262147 7.230256 7.170210

[345] 7.137925 7.104712 7.121480 7.118843 7.117378 7.090968 7.084255 7.108023

[353] 7.051757 7.028121 7.054859 7.078414 7.114706 7.079548 7.114003 7.078881

[361] 7.090100 7.129850 7.126992 7.139105 7.094334 7.087728 7.105869 7.068680

[369] 7.077331 7.074386 7.037065 7.015845 7.031067 7.045745 7.045352 7.067747

[377] 7.079858 7.077529 7.085807 7.101523 7.146237 7.158741 7.093481 7.064842

[385] 7.096556 7.130847 7.141951 7.115001 7.112850 7.150857 7.145183 7.094028

[393] 7.118288 7.150233 7.145857 7.144554 7.146988 7.146072 7.176767 7.107342

[401] 7.123751 7.149823 7.156241 7.111539 7.146516 7.174778 7.180236 7.151633

[409] 7.177896 7.097421 7.134737 7.139064 7.180692 7.176342 7.181836 7.204021

[417] 7.208299 7.146985 7.148830 7.107342 7.088155 7.052253 7.006641 7.015996

[425] 7.046809 7.061120 7.032716 7.041052 7.064903 7.059522 7.091234 7.090501

[433] 7.081098 7.045878 7.059936 7.085968 7.130180 7.102629 7.103868 7.103656

[441] 7.117224 7.131319 7.098237 7.118564 7.101916 7.079309 7.078102 7.027398

[449] 7.017138 7.023637 7.091488 7.122048 7.133928 7.131692 7.159420 7.169551

[457] 7.146890 7.142987 7.140642 7.126635 7.175829 7.154942 7.128365 7.181859

[465] 7.175072 7.169410 7.238363 7.237141 7.289526 7.306859 7.273925 7.263208

[473] 7.243585 7.268352 7.224143 7.250465 7.270233 7.255779 7.240344 7.291077

[481] 7.307398 7.266715 7.266844 7.284025 7.303616 7.278101 7.329640 7.312715

[489] 7.345105 7.305630 7.276796 7.269477 7.256781 7.295397 7.296595 7.320085

[497] 7.331981 7.331155 7.333449 7.307147 7.316371 7.328005 7.300980 7.281100

[505] 7.288726 7.225282 7.232577 7.229973 7.249190 7.323170 7.325466 7.359743

[513] 7.354087 7.314519 7.316979 7.316903 7.295794 7.266593 7.290003 7.309931

[521] 7.269465 7.240654 7.244580 7.236161 7.220369 7.216624 7.214971 7.182581

[529] 7.164027 7.127042 7.135585 7.153463 7.168115 7.128040 7.155798 7.135208

[537] 7.177279 7.221284 7.256975 7.279565 7.272557 7.250464 7.257887 7.278663

[545] 7.233752 7.227740 7.232923 7.291705 7.251748 7.216629 7.246354 7.246186

[553] 7.222686 7.240670 7.240207 7.191035 7.232885 7.270038 7.227842 7.217455

[561] 7.170541 7.178615 7.206536 7.192185 7.185152 7.186717 7.157135 7.174848

[569] 7.183865 7.156051 7.142732 7.105314 7.068098 7.087493 7.020719 7.031725

[577] 7.016702 7.015673 6.999963 7.045578 7.069696 7.110356 7.116568 7.111650

[585] 7.046876 7.098489 7.098522 7.104072 7.135575 7.074844 7.090479 7.100506

[593] 7.121507 7.123277 7.090857 7.085827 7.079468 7.036278 7.043417 7.009129

[601] 7.006092 6.995034 6.997751 6.961613 6.962277 7.000038 6.992817 7.028478

[609] 7.030588 7.004229 6.975593 6.957772 6.948820 6.938782 6.912071 6.900986

[617] 6.913416 6.942149 6.958355 7.004749 6.981839 7.010009 7.034225 7.090857

[625] 7.087262 7.141156 7.165744 7.170199 7.123199 7.155267 7.092919 7.153585

[633] 7.159173 7.095379 7.123217 7.091858 7.103740 7.071565 7.035169 7.066140

[641] 7.095298 7.084756 7.093266 7.078597 7.089700 7.139309 7.117069 7.137335

[649] 7.136439 7.116762 7.173263 7.133395 7.147188 7.203146 7.242657 7.250927

[657] 7.201513 7.188890 7.169367 7.157233 7.103280 7.040244 7.070839 7.083728

[665] 7.058956 7.024028 7.031458 7.069782 7.106274 7.130350 7.126406 7.142285

[673] 7.132842 7.136819 7.126732 7.165840 7.160837 7.147163 7.121239 7.107829

[681] 7.094635 7.127397 7.129823 7.144622 7.129571 7.127359 7.095176 7.072701

[689] 7.025348 7.026530 7.053881 7.038662 7.077207 7.065855 7.092630 7.090197

[697] 7.104571 7.065959 7.086947 7.086198 7.061146 7.071278 7.087874 7.075494

[705] 7.084810 7.071096 7.017197 7.028904 7.045269 6.995492 6.993004 7.002403

[713] 7.070153 7.048646 7.053306 7.068167 7.065019 7.051720 7.029827 7.067742

[721] 7.101737 7.074833 7.109695 7.117952 7.171604 7.157373 7.122277 7.102073

[729] 7.122280 7.092120 7.065516 7.036670 7.006340 7.011977 7.011673 6.996891

[737] 7.014798 7.015064 6.995969 6.958202 6.973073 6.979068 6.937016 6.969629

[745] 6.984088 6.976605 6.951565 6.942892 6.924040 6.907088 6.912183 6.905260

[753] 6.903310 6.956582 6.930251 6.886287 6.911363 6.924639 6.914682 6.911524

[761] 6.893737 6.899142 6.907357 6.914254 6.927394 6.909583 6.892624 6.903670

[769] 6.903444 6.902892 6.918966 6.955622 6.950996 6.933009 6.895144 6.931373

[777] 6.889020 6.882409 6.874159 6.888222 6.897185 6.883539 6.915821 6.932445

[785] 6.940479 6.929080 6.931299 6.964280 6.999782 6.989388 6.981434 6.968265

[793] 7.003788 6.995484 7.003156 6.954861 6.982269 6.972628 6.998854 7.000844

[801] 6.959610 6.940751 6.952307 6.919167 6.934823 6.991495 7.009198 7.024987

[809] 7.010942 7.009779 7.026069 7.053679 7.069666 7.098924 7.100481 7.088994

[817] 7.088969 7.037677 7.023283 7.052647 7.074681 7.051656 7.132211 7.057101

[825] 7.065356 7.075734 7.102188 7.155825 7.167115 7.153556 7.088877 7.113851

[833] 7.099440 7.106808 7.036434 7.001959 7.003078 7.086281 7.089996 7.053244

[841] 7.093258

> theta.accept

[1] 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 0 0 0 1 1 1 0 1 1 1 1

[37] 0 1 0 1 1 1 0 0 1 1 1 1 0 0 1 0 1 1 0 1 0 1 0 1 1 0 1 0 1 1 1 1 0 1 1 1

[73] 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1

[109] 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1

[145] 0 1 1 0 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1

[181] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1

[217] 1 1 1 1 1 0 1 1 0 0 0 0 1 1 1 0 1 1 0 0 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1

[253] 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[289] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 1 1

[325] 1 0 1 1 1 1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1

[361] 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1

[397] 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 0

[433] 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 0

[469] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1

[505] 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1

[541] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1

[577] 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 1 1 0

[613] 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[649] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[685] 1 1 1 0 0 1 0 1 1 1 0 1 1 1 0 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1

[721] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1

[757] 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 1 1 1 1 1

[793] 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 0 0 1 0 0 1 1 0 1 1 1 0 1 1

[829] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1

[865] 0 1 1 1 1 1 1 0 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1

[901] 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 1 1 0

[937] 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 1 1 1

[973] 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1

[ reached getOption("max.print") -- omitted 500 entries ]

> colchain <- rep("blue",nchain)

> for (ichain in 1: nchain){

+ if (theta.accept[ichain] == 0) colchain[ichain]="red"

+ }

> # Traceplots

> # Figure 6.13

> par(mar=c(4,5,4,2))

> par(mfrow=c(1,1))

> par(pty="m")

> plot(nstart:nchain,theta.1[nstart:nchain],xlab="Iteration",

+ ylab=expression(paste(mu)),cex.lab=1.5,cex=1.3,cex.axis=1.3,

+ type="n",main="",bty="o")

> for (ichain in nstart:(nchain-1)){

+ segments(ichain,theta.1[ichain],ichain+1,theta.1[ichain+1],

+ col=colchain[ichain+1],lwd=2)

+ }

> text(500,7.3,"(a)",cex=1.5)

> plot(nstart:nchain,theta.2[nstart:nchain],xlab="Iteration",

+ ylab=expression(paste(sigma^{2})),cex.lab=1.5,cex=1.5,cex.axis=1.3,

+ type="n",main="",bty="o")

> for (ichain in nstart:(nchain-1)){

+ segments(ichain,theta.2[ichain],ichain+1,theta.2[ichain+1],

+ col=colchain[ichain+1],lwd=2)

+ }

> text(500,2.37,"(b)",cex=1.5)

> # Marginal posterior distributions of mu and signa^2
> # After burn-in part
> 
> unl.mu <- 6.8
> upl.mu <- 7.4
> 
> unl.sigma2 <- 1.2
> upl.sigma2 <- 2.5
> 
> 
> nmu <- 50
> nsigma2 <- 50
> ntilde <- 50
> 
> # Construction of true posteriors
> 
> mu <- seq(unl.mu,upl.mu,length=nmu)
> mus <- (mu-mtalp)/(sdtalp/sqrt(ntalp))
> pmu <- dt(mus,ntalp-1)*sqrt(ntalp)/sdtalp
> sigma2 <- seq(unl.sigma2,upl.sigma2,length=nsigma2)
> sigma2s <- (ntalp-1)*vartalp/sigma2
> psigma2 <- dchisq(sigma2s,ntalp-1)*(ntalp-1)*vartalp/(sigma2^2)
> 
> # Show sampled posteriors obtained from MH algorithm
> 
> par(mfrow=c(1,2))
> par(pty="s")
> hist(theta.1[nstart:nchain],xlab=expression(mu),ylab="", 
+      nclas=30,prob=T,main="",cex.lab=1.5,cex.axis=1.3,col="light blue")
> lines(mu,pmu,lwd=3,lty=1,col="black")
> text(6.9,6,"(a)",cex=1.3)
> 
> hist(theta.2[nstart:nchain],xlab=expression(paste(sigma^{2})),ylab="", 
+      nclas=30,prob=T,main="",cex.lab=1.5,cex.axis=1.3,col="light blue")
> lines(sigma2,psigma2,lwd=3,lty=1,col="black")
> text(1.55,3.5,"(b)",cex=1.3)
> 
> summary(theta.1[nstart:nchain]);summary(theta.2[nstart:nchain])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.874   7.052   7.110   7.114   7.170   7.360 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.525   1.731   1.819   1.822   1.896   2.204 
>

> # Marginal posterior distributions of mu and signa^2

> # After burn-in part

> unl.mu <- 6.8

> upl.mu <- 7.4

> unl.sigma2 <- 1.2

> upl.sigma2 <- 2.5

> nmu <- 50

> nsigma2 <- 50

> ntilde <- 50

> # Construction of true posteriors

> mu <- seq(unl.mu,upl.mu,length=nmu)

> mus <- (mu-mtalp)/(sdtalp/sqrt(ntalp))

> pmu <- dt(mus,ntalp-1)*sqrt(ntalp)/sdtalp

> sigma2 <- seq(unl.sigma2,upl.sigma2,length=nsigma2)

> sigma2s <- (ntalp-1)*vartalp/sigma2

> psigma2 <- dchisq(sigma2s,ntalp-1)*(ntalp-1)*vartalp/(sigma2^2)

> # Show sampled posteriors obtained from MH algorithm

> par(mfrow=c(1,2))

> par(pty="s")

> hist(theta.1[nstart:nchain],xlab=expression(mu),ylab="",

+ nclas=30,prob=T,main="",cex.lab=1.5,cex.axis=1.3,col="light blue")

> lines(mu,pmu,lwd=3,lty=1,col="black")

> text(6.9,6,"(a)",cex=1.3)

> hist(theta.2[nstart:nchain],xlab=expression(paste(sigma^{2})),ylab="",

+ nclas=30,prob=T,main="",cex.lab=1.5,cex.axis=1.3,col="light blue")

> lines(sigma2,psigma2,lwd=3,lty=1,col="black")

> text(1.55,3.5,"(b)",cex=1.3)

> summary(theta.1[nstart:nchain]);summary(theta.2[nstart:nchain])

Min. 1st Qu. Median Mean 3rd Qu. Max.

6.874 7.052 7.110 7.114 7.170 7.360

Min. 1st Qu. Median Mean 3rd Qu. Max.

1.525 1.731 1.819 1.822 1.896 2.204

Science To Medicine

Just My Daily Study Note by ts.anesth.kpum

メトロポリス・ヘイスティング・アルゴリズム