[Lecture Notes]Topics in Bayesian statistics

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

本帖隐藏的内容

Topics in Bayesian statistics.pdf (931.21 KB)

2016-6-5 01:42:33 上传

http://people.bath.ac.uk/masss/ma40189.html

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：Statistics statistic Bayesian Lecture Statist

1. gibbs.update=function(start=c(1,1), n=10, rho=0.5, pause=T) {
   
2. # use pause=T for small n, to display step by step pressing return
   
3. # in between.
   
4. # use pause-F to run stragiht through and plot the pictures only
   
5. # at the end
   
6. 
   
7. # set up vectors for the samples
   
8.   x=rep(NA,n)
   
9.   y=rep(NA,n)
   
10.   x[1]=start[1]
    
11.   y[1]=start[2]
    
12. 
    
13.   #Set up plots
    
14.   xx=seq(-3,3,len=20)
    
15.   yy=seq(-3,3,len=20)
    
16.   z = outer(xx,yy,f,rho)
    
17.   if (pause) {
    
18.     layout(matrix(c(0,1,1,0,0,1,1,0,2,2,3,3,2,2,3,3),4,4,byrow=T))
    
19.     #plot contours
    
20.     contour(xx,yy,z,col="blue",levels=c(0.01,0.02,0.05,0.1,0.14),
    
21.             drawlabels=F)
    
22.     title(xlab="x",ylab="y")
    
23.     #mark starting point:
    
24.     points(x[1],y[1],pch=4,col="red")
    
25.     #Time series plot so far:
    
26.     # X:
    
27.     plot(c(0,n),c(min(start[1],-3),max(start[1],3)),
    
28.          xlab="iterations",ylab="x",type="n")
    
29.     lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
30.     lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
31.     # Y:
    
32.     plot(c(0,n),c(min(start[2],-3),max(start[2],3)),
    
33.          xlab="iterations",ylab="y",type="n")
    
34.     lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
35.     lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
36.     # Wait for key press
    
37.     readLines(n=1)
    
38.   }
    
39. 
    
40.    
    
41.   #Main loop
    
42.   for (t in 2:n) {
    
43. 
    
44.     # sample x from p(x|y=y[t-1]):
    
45.     x[t]=rnorm(1, rho* y[t-1], sqrt(1-rho^2))
    
46.     # sample y from p(y|=x[t]):
    
47.     y[t]=rnorm(1, rho* x[t], sqrt(1-rho^2))
    
48. 
    
49.     if (pause || (t==n)) {
    
50.    
    
51.       #plot contours (use xx,yy and z from before)
    
52.       contour(xx,yy,z ,col="blue",levels=c(0.01,0.02,0.05,0.1,0.14),
    
53.               drawlabels=F)
    
54.       title(xlab="x",ylab="y")
    
55.       #mark starting point:
    
56.       points(x[1],y[1],pch=4,col="red")
    
57.       #draw path of the chain
    
58.       lines(c(x[1],rep(x[2:t],rep(2,t-1))),
    
59.             c(rep(y[1:(t-1)],rep(2,t-1)),y[t]),col="red")
    
60. 
    
61.       #Time series plot so far:
    
62.       # X:
    
63.       plot(c(0,n),c(min(start[1],-3),max(start[1],3)),
    
64.            xlab="iterations",ylab="x",type="n")
    
65.       lines(1:t,x[1:t],col="red")
    
66.       lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
67.       lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
68.       # Y:
    
69.       plot(c(0,n),c(min(start[2],-3),max(start[2],3)),
    
70.            xlab="iterations",ylab="y",type="n")
    
71.       lines(1:t,y[1:t],col="red")
    
72.       lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
73.       lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
74.       
    
75.     }
    
76. 
    
77.     if ((pause)  && (t<n)) {
    
78.       # Wait for key press
    
79.       readLines(n=1)
    
80.     }
    
81. 
    
82. 
    
83.   #end main loop
    
84.   }
    
85. }

复制代码

1. metropolis = function (start=1, n=100, sig.q=1, mu.p=0, sig.p=1) {
   
2. 
   
3. #starting value theta_0 = start
   
4. #run for n iterations
   
5. 
   
6.   #set up vector theta[] for the samples
   
7.   theta=rep(NA,n)
   
8.   theta[1]=start
   
9. 
   
10.   #count the number of accepted proposals
    
11.   accept=0
    
12. 
    
13.   #Main loop
    
14.   for (t in 2:n) {
    
15. 
    
16.     #pick a candidate value theta* from N(theta[t-1], sig.q)
    
17.     theta.star = rnorm(1, theta[t-1], sqrt(sig.q))
    
18.    
    
19.     #generate a random number between 0 and 1:
    
20.     u = runif(1)
    
21. 
    
22.     # calculate acceptance probability:
    
23.     r = (dnorm(theta.star,mu.p,sig.p)) / (dnorm(theta[t-1],mu.p,sig.p) )
    
24.    
    
25.     a=min(r,1)
    
26.    
    
27. 
    
28.     #if u<a we accept the proposal; otherwise we reject
    
29.     if (u<a) {   
    
30.         #accept
    
31.             accept=accept+1
    
32.             theta[t] = theta.star
    
33.     }
    
34.     else {
    
35.         #reject
    
36.         theta[t] = theta[t-1]
    
37.         }
    
38.   #end loop
    
39.   }

复制代码

1. metropolis.update =function(start=1, n=10, sig.q=1, mu.p=0, sig.p=1) {
   
2. 
   
3.   #set up vector theta[] for the samples
   
4.   theta=rep(NA,n)
   
5.   theta[1]=start
   
6. 
   
7.   #count the number of accepted proposals
   
8.   accept=0
   
9. 
   
10.   #Set up plots
    
11.   layout(matrix(c(0,1,1,0,0,1,1,0,2,2,2,2,2,2,2,2),4,4,byrow=T))
    
12.   #Selection of points
    
13.   plot (c(-5,5),c(0,0.5),type="n",main="",xlab="theta",ylab="density")
    
14.   #Target density
    
15.   curve(dnorm(x,mu.p,sqrt(sig.p)),add=T,lty=2,col="blue")
    
16.   #Proposal density
    
17.   curve(dnorm(x,theta[1],sqrt(sig.q)),add=T,col="red")
    
18.   #current proposal
    
19.   lines(c(theta[1],theta[1]),c(0,0.5),col="red")
    
20. 
    
21.   #Time series plot so far:
    
22.   plot(c(0,n),c(min(start,-3),max(start,3)),
    
23.        xlab="iterations",ylab="theta",type="n")
    
24.   lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
25.   lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
26. 
    
27.   readLines(n=1)
    
28.   
    
29.   #Main loop
    
30.   for (t in 2:n) {
    
31. 
    
32.     #pick a candidate value theta* from N(theta[t-1], sig.q)
    
33.     theta.star = rnorm(1, theta[t-1], sqrt(sig.q))
    
34.    
    
35.     #generate a random number between 0 and 1:
    
36.     u = runif(1)
    
37. 
    
38.     # calculate acceptance probability:
    
39.     r = (dnorm(theta.star,mu.p,sig.p)) / (dnorm(theta[t-1],mu.p,sig.p) )
    
40.     a=min(r,1)
    
41.    
    
42. 
    
43.     #if u<a we accept the proposal; otherwise we reject
    
44.     if (u<a) {   
    
45.         #accept
    
46.             accept=accept+1
    
47.             theta[t] = theta.star
    
48.             cat("\ttheta*=",theta.star,"\n")
    
49.             cat("ACCEPT: theta[t]=",theta[t],"\n")
    
50.     }
    
51.     else {
    
52.         #reject
    
53.         theta[t] = theta[t-1]
    
54.         cat("\ttheta*=",theta.star,"\n")
    
55.         cat("REJECT: theta[t]=",theta[t],"\n")
    
56.         }
    
57. 
    
58.     # update plots
    
59.     layout(matrix(c(0,1,1,0,0,1,1,0,2,2,2,2,2,2,2,2),4,4,byrow=T))
    
60.     #Selection of points
    
61.     plot (c(-5,5),c(0,0.5),type="n",main="",xlab="theta",ylab="density")
    
62.     #Target density
    
63.     curve(dnorm(x,mu.p,sqrt(sig.p)),add=T,lty=2,col="blue")
    
64.     #Proposal density
    
65.     curve(dnorm(x,theta[t-1],sqrt(sig.q)),add=T,col="red")
    
66.     #current proposal
    
67.     lines(c(theta.star,theta.star),c(0,0.5),col="red",lty=2)
    
68.     readLines(n=1)
    
69.     #theta[t]
    
70.     lines(c(theta[t],theta[t]),c(0,0.5),col="red")
    
71. 
    
72.     #Plot the density so far
    
73.     if (t>2) lines(density(theta[1:t]),col="green",lty=2)
    
74.     #Time series plot so far:
    
75.     plot(c(0,n),c(min(start,-3),max(start,3)),
    
76.          xlab="iterations",ylab="theta",type="n")
    
77.     lines(0:t,theta[1:(t+1)],col="red")
    
78.     lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    
79.     lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    
80. 
    
81.       
    
82.     #Wait for keypress
    
83.     readLines(n=1)
    
84.    
    
85.   #end loop
    
86.   }
    
87. 
    
88. 
    
89. }

复制代码

Metropolis-Hastings for Bivariate Normals

1. metropolis1=function( n=1000 , rho=0.99 , start1=-1 , start2=-1 , tau=1 )
   
2. {
   
3.   x1 <- c(1:n)
   
4.   x2 <- c(1:n)
   
5.   x1[1] <- start1
   
6.   x2[1] <- start2
   
7.   par( mfrow=c( 2,2 ) )
   
8.   for ( i in 2:n )
   
9.    {
   
10.     try <- x1[i-1] + rnorm( 1 , 0 , tau )
    
11.     temp <- exp( ( try^2 - x1[i-1]^2 - 2.0 * rho * x2[i-1] * ( try - x1[i-1] ) )/(-2.0*(1-rho^2)) )
    
12.     if ( runif(1) < temp )
    
13.      {
    
14.       x1[i] <- try
    
15.      }
    
16.     else
    
17.      {
    
18.       x1[i] <- x1[i-1]
    
19.      }
    
20.     try <- x2[i-1] + rnorm(1,0,tau )
    
21.     temp <- exp( ( try^2 - x2[i-1]^2 - 2.0 * rho * x1[i] *( try - x2[i-1]) )/(-2.0*(1-rho^2)) )
    
22.     if ( runif(1) < temp )
    
23.      {
    
24.       x2[i] <- try
    
25.      }
    
26.     else
    
27.      {
    
28.       x2[i] <- x2[i-1]
    
29.      }
    
30. 
    
31. 
    
32.   plot( x1[1:i] , type="l" , xlab="Iteration" ,xlim=c(1,n),ylab="x1")
    
33.   plot( x2[1:i] , type="l" , xlab="Iteration",xlim=c(1,n),ylab="x2" )
    
34.   plot( x1[1:i] , x2[1:i] , type="p" , pch="." ,xlab="x1",ylab="x2")
    
35.   plot( x1[1:i] , x2[1:i] , type="l" ,col="red",xlab="x1",ylab="x2")
    
36. }
    
37. }

复制代码

Metropolis-Hastings for sampling from a mixture of bivariate normals

1. metropolis2=function( n=1000 , rho=0.0 , start1=0 , start2=0 , tau=0.25 , sigma2=0.1 )
   
2. {
   
3.   x1 <- c(1:n)
   
4.   x2 <- c(1:n)
   
5.   x1[1] <- start1
   
6.   x2[1] <- start2
   
7.   scale <- -0.5 / ( sigma2* ( 1 - rho^2 ) )
   
8.   breaks <- c( -20:(1+20) )/4
   
9.   par( mfrow=c( 2,2 ) )
   
10.   for ( i in 2:n )
    
11.    {
    
12.     try <- x1[i-1] + rnorm( 1 , 0 , tau )
    
13.     temp <- exp( scale * ( try^2 - 2*rho*try*x2[i-1] +x2[i-1]^2) ) +
    
14.             exp( scale * ( (try-1)^2 - 2*rho*(try-1)*(x2[i-1]-1) + (x2[i-1]-1)^2 ) )
    
15.     temp <- temp / ( exp( scale * ( x1[i-1]^2 - 2*rho*x1[i-1]*x2[i-1] + x2[i-1]^2 ) ) +
    
16.                      exp( scale * ( (x1[i-1]-1)^2 - 2*rho*(x1[i-1]-1)*(x2[i-1]-1) + (x2[i-1]-1)^2 ) ) )
    
17.     if ( runif(1) < temp )
    
18.      {
    
19.       x1[i] <- try
    
20.      }
    
21.     else
    
22.      {
    
23.       x1[i] <- x1[i-1]
    
24.      }
    
25.     try <- x2[i-1] + rnorm(1,0,tau )
    
26.     temp <- exp( scale * ( try^2 - 2*rho*try*x1[i] + x1[i]^2) ) +
    
27.             exp( scale * ( (try-1)^2 - 2*rho*(try-1)*(x1[i]-1) + (x1[i]-1)^2 ) )
    
28.     temp <- temp / ( exp( scale * ( x2[i-1]^2 - 2*rho*x1[i]*x2[i-1] + x1[i]^2 ) ) +
    
29.                      exp( scale * ( (x2[i-1]-1)^2 - 2*rho*(x1[i]-1)*(x2[i-1]-1) +(x1[i]-1)^2  ) ) )
    
30.     if ( runif(1) < temp )
    
31.      {
    
32.       x2[i] <- try
    
33.      }
    
34.     else
    
35.      {
    
36.       x2[i] <- x2[i-1]
    
37.      }
    
38. 
    
39. 
    
40.   hist( x1[1:i] , breaks , main="x1" , xlab="x1")
    
41.   hist( x2[1:i] , breaks , main="x2" , xlab="x2" )
    
42.   plot( x1[1:i] , x2[1:i] , type="p" , pch="." ,xlab="x1",ylab="x2")
    
43.   plot( x1[1:i] , x2[1:i] , type="l" ,col="red",xlab="x1",ylab="x2")
    
44. }
    
45. }

复制代码

Gibbs Sampler from Bivariate Normals

1. gibbs1=function( n=1000 , rho=0.99 , start1=-1 , start2=-1 )
   
2. {
   
3.   x1 <- c(1:n)
   
4.   x2 <- c(1:n)
   
5.   x1[1] <- start1
   
6.   x2[1] <- start2
   
7.   par( mfrow=c( 2,2 ) )
   
8.   for ( i in 2:n )
   
9.    {
   
10.     x1[i] <- rnorm( 1 , x2[i-1] * rho , 1-rho^2 )
    
11.     x2[i] <- rnorm( 1, x1[i] * rho , 1 - rho^2 )
    
12. 
    
13.   plot( x1[1:i] , type="l" , xlab="Iteration" ,xlim=c(1,n),ylab="x1")
    
14.   plot( x2[1:i] , type="l" , xlab="Iteration",xlim=c(1,n),ylab="x2" )
    
15.   plot( x1[1:i] , x2[1:i] , type="p" , pch="." ,xlab="x1",ylab="x2")
    
16.   plot( x1[1:i] , x2[1:i] , type="l" ,col="red",xlab="x1",ylab="x2")
    
17. }
    
18. }

复制代码

谢谢分享

感恩分享