1ordergibbs=function(data,m)
2{
3# implements Gibbs sampling for table of means
4# with prior belief in order restriction
5# input:  data = data matrix with two columns [sample mean, sample size]
6#         m = number of iterations of Gibbs sampling
7# output:  matrix of simulated values of means where each row
8#          represents one simulated draw
9
10#####################################################
11rnormt=function(n,mu,sigma,lo,hi)
12{
13# simulates n random variates from a normal(mu,sigma)
14# distribution truncated on the interval (lo, hi)
15
16p=pnorm(c(lo,hi),mu,sigma)
17return(mu+sigma*qnorm(runif(n)*(p[2]-p[1])+p[1]))
18}
19#####################################################
20
21y=data[,1]  # sample means
22n=data[,2]  # sample sizes
23s=.65       # assumed value of sigma for this example
24I=8; J=5    # number of rows and columns in matrix
25
26# placing vectors y, n into matrices
27
28y=t(array(y,c(J,I)))
29n=t(array(n,c(J,I)))
30y=y[seq(8,1,by=-1),]
31n=n[seq(8,1,by=-1),]
32
33# setting up the matrix of values of the population means mu
34# two rows and two columns are added that help in the simulation
35# of individual values of mu from truncated normal distributions
36
37mu0=Inf*array(1,c(I+2,J+2))
38mu0[1,]=-mu0[1,]
39mu0[,1]=-mu0[,1]
40mu0[1,1]=-mu0[1,1]
41mu=mu0
42
43# starting value of mu that satisfies order restriction
44
45m1=c(2.64,3.02,3.02,3.07,3.34)
46m2=c(2.37,2.63,2.74,2.76,2.91)
47m3=c(2.37,2.47,2.64,2.66,2.66)
48m4=c(2.31,2.33,2.33,2.33,2.33)
49m5=c(2.04,2.11,2.11,2.33,2.33)
50m6=c(1.85,1.85,1.85,2.10,2.10)
51m7=c(1.85,1.85,1.85,1.88,1.88)
52m8=c(1.59,1.59,1.59,1.67,1.88)
53muint=rbind(m8,m7,m6,m5,m4,m3,m2,m1)
54mu[2:(I+1),2:(J+1)]=muint
55
56MU=array(0,c(m,I*J))  # arry MU stores simulated values of mu
57
58#####################  main loop #######################
59for (k in 1:m)
60{
61	for (i in 2:(I+1))
62	{
63		for (j in 2:(J+1))
64		{
65			lo=max(c(mu[i-1,j],mu[i,j-1]))
66			hi=min(c(mu[i+1,j],mu[i,j+1]))
67			mu[i,j]=rnormt(1,y[i-1,j-1],s/sqrt(n[i-1,j-1]),lo,hi)
68		}
69	}
70	mm=mu[2:(I+1),2:(J+1)]
71	MU[k,]=array(mm,c(1,I*J))
72}
73
74return(MU)
75}
76