# Uploaded: April 09, 2007

#######################

avgprob <- function (n, r, measure, samp.prior, param1, param2, 
    delta = 0.80, tdi.p0 = 0.80, alpha = 0.05, AB = 0.001, 
    n.rep = 100, n.iter = 1000, n.burn = 250)
{
    if (identical(measure, "ccc")) {
        t.nr <- replicate(n.rep, ccc.prob (n, r, samp.prior, param1, param2,
        delta, alpha, AB, n.iter, n.burn))
        res <- apply(t.nr, 1, mean)
        names(res) <- c("t11star", "t22star", "t12star")
        return(res)
    }
    else if  (identical(measure, "tdi")) {
        t.nr <- replicate(n.rep, tdi.prob (n, r, samp.prior, param1, param2, 
        delta, tdi.p0, alpha, AB, n.iter, n.burn))
        res <- apply(t.nr, 1, mean)
        names(res) <- c("t11star", "t22star", "t12star")
        return(res)
    }
    else (stop (" 'measure' must be either 'ccc' or 'tdi' "))
    
}

####################

ccc.prob <- function(n, r, samp.prior, param1, param2,
        delta, alpha, AB, n.iter, n.burn)
{

# Computes posterior probabilities for concordance correlation
# as the measure of agreement

    if ( identical(samp.prior, "normal")) {
        mu <- rnorm(2, mean=param1[1:2], sd=param2[1:2])
        theta <- exp(rnorm(4, mean=param1[3:6], sd=param2[3:6]))
    }
    else if ( identical(samp.prior, "uniform")) {
        mu <- runif(2, min=param1[1:2], max=param2[1:2])
        theta <- exp(runif(4, min=param1[3:6], max=param2[3:6]))                       
    }
    else (stop (" 'samp.prior' must be either 'uniform' or 'normal' "))
    
    if( identical(r, 1)) theta[4] <- 0

    Y <- simData(mu, theta, n, r)
    
    sim.post <- simPost(AB, theta, Y, n.iter, n.burn)
    mu.post <- sim.post$mu.post
    theta.post <- sim.post$theta.post

    ccc.post <- ccc.fun(mu.post[,1], mu.post[,2], 
                    theta.post[,1], theta.post[,2], 
                    theta.post[,3], theta.post[,4])
                    
    cc11 <- ccc.post$intra1
    cc22 <- ccc.post$intra2
    cc12 <- ccc.post$inter

    t11 <- -log(cc11)
    t22 <- -log(cc22)
    t12 <- -log(cc12)

    t.ubd <- apply(cbind(t11, t22, t12),2,quantile,prob=(1-alpha),type=4)
    ubd.delta<-delta*t.ubd

    F11 <- ecdf(t11)
    F22 <- ecdf(t22)
    F12 <- ecdf(t12)
    
    prob <- c(F11(ubd.delta[1]), F22(ubd.delta[2]), F12(ubd.delta[3]))   
    return(prob)
} 

##

ccc.fun <- function(mu1, mu2, theta1, theta2, theta3, theta4) 
{

var1 <- theta3 + theta4 + theta1
var2 <- theta3 + theta4 + theta2

intra1 <- (theta3 + theta4)/var1
intra2 <- (theta3 + theta4)/var2
inter <- (2*theta3)/( var1 + var2 + (mu1-mu2)^2)

return(list(intra1=intra1, intra2=intra2, inter=inter))
}

####################

tdi.prob <- function(n, r, samp.prior, param1, param2, delta, tdi.p0, 
    alpha, AB, n.iter, n.burn) 
{

# Computes posterior probabilities for total deviation index
# as the measure of agreement


# mu = (beta1, beta2)
# theta = (sigma1^2, sigma2^2, psi^2, psi_I^2)

    if ( identical(samp.prior, "normal")) {
        mu <- rnorm(2, mean=param1[1:2], sd=param2[1:2])
        theta <- exp(rnorm(4, mean=param1[3:6], sd=param2[3:6]))
    }
    else if ( identical(samp.prior, "uniform")) {
        mu <- runif(2, min=param1[1:2], max=param2[1:2])
        theta <- exp(runif(4, min=param1[3:6], max=param2[3:6]))                       
    }
    else (stop (" 'samp.prior' must be either 'uniform' or 'normal' "))
    
    if( identical(r, 1)) theta[4] <- 0

    Y <- simData(mu, theta, n, r)
    
    sim.post <- simPost(AB, theta, Y, n.iter, n.burn)          
    mu.post <- sim.post$mu.post
    theta.post <- sim.post$theta.post
    
    mu.p <- mu.post[,1]-mu.post[,2]
    sigma.p <- sqrt(theta.post[,1]+theta.post[,2]+2*theta.post[,4])

    q11 <- sqrt(theta.post[,1])*sqrt(2*qchisq(tdi.p0,1))
    q22 <- sqrt(theta.post[,2])*sqrt(2*qchisq(tdi.p0,1))
    q12 <- q.p(tdi.p0, mu.p, sigma.p)

    ubd.q <- apply(cbind(q11, q22, q12),2,quantile,prob=(1-alpha),type=4)
    ubd.delta<-delta*ubd.q
 
    F11 <- ecdf(q11)
    F22 <- ecdf(q22)
    F12 <- ecdf(q12)
 
    prob <- c(F11(ubd.delta[1]), F22(ubd.delta[2]), F12(ubd.delta[3]))      
    return(prob)
} 

##

q.p <- function(p, du,sdev) 
{
#calculates the p_th qunatile of |diff|
#du and sdev can be vector, but must have the same length
    ans <- mapply(q.scalar, m=du, s=sdev, p=p)
    return(ans)
}

##

q.scalar <- function(p, m, s) 
{
# m=mean,s=sigma, p=prob
#p must in (0.5,1)

    a<-m/s
    p<-1+p
    llim <- abs(a)
    ulim <- llim + 1.e-6 + qnorm(p/2)

    f <- function(q) { pnorm(q-a)+pnorm(q+a)- p}
    
    ans <- uniroot(f, c(llim, ulim), tol=1.e-6, maxiter=1e+4)$root
    return(s*ans)
}

#############################

simData <- function(mu,theta,n,r) 
{
#simulates data for given parameter values.
    sigma<-sqrt(theta)
    if(identical(r, 1)) sigma[4] <- 0
    bI <- rnorm(n)*sigma[3]
    Y1 <- mu[1]+bI+rnorm(n)*sigma[4]+rnorm(n*r)*sigma[1]
    Y2 <- mu[2]+bI+rnorm(n)*sigma[4]+rnorm(n*r)*sigma[2]
    return(list(Y1=matrix(Y1,ncol=r),Y2=matrix(Y2,ncol=r)))
}


################################
   
simPost <- function(AB, initTheta, data, n.iter, n.burn)
{
#simulates from posterior using Gibbs sampler.
         
    theta<-initTheta      
    stat<-lapply(data,data.stat)
    theta.par<-list(rMeans=rbind(stat$Y1$rMeans,stat$Y2$rMeans),
        sv=c(stat$Y1$sv,stat$Y2$sv),dim=stat$Y1$dim)
    dY<-rbind(stat$Y1$rMeans.c, stat$Y2$rMeans.c, 
        (stat$Y1$rMeans.c-stat$Y2$rMeans.c))
    gamma.par<-list(mean=c(stat$Y1$mean,stat$Y2$mean),dY=dY,dim=stat$Y1$dim)
        
    for(burn in seq(n.burn)){
         gamma<-gamma.p(theta,gamma.par)
         theta<-theta.p(gamma,theta.par,AB)
    }
    post.draw<-n.iter-n.burn
    q.emf<-seq(post.draw)
    mu.post<-matrix(0,nrow=post.draw,ncol=2)
    theta.post<-matrix(0,nrow=post.draw,ncol=4)
   for( size in seq(post.draw)){
      gamma<-gamma.p(theta,gamma.par)
      theta<-theta.p(gamma,theta.par,AB)
      mu.post[size,]<-gamma$mu
      theta.post[size,]<-theta
     }         
     return(list(mu.post=mu.post,theta.post=theta.post))
}

#################################

data.stat <- function(data)
{
# computes summary statistics for the data
   Y1<-data
   Y1.dim<-dim(Y1)  #n*r  
   Y1.rm<-rowMeans(Y1)
   Y1.sv<-sum((Y1-Y1.rm)^2)
   Y1.m<-mean(Y1.rm)
   Y1.c<-Y1.rm-Y1.m
   return(list(mean=Y1.m,sv=Y1.sv,dim=Y1.dim,rMeans=Y1.rm,rMeans.c=Y1.c))
}

###############################

theta.p <- function(gamma,parY,AB)
{
#computes full conditional of theta given gamma(=(mu,b)) and csY
    A<-rep(AB,4)
    B<-A
    
    n<-parY$dim[1] 
    r<-parY$dim[2]
    
    mu<-gamma$mu
    bs<-gamma$bs
    
    dYb1<-parY$rMeans[1,]-mu[1]-bs[1,]-bs[2,]
    dYb2<-parY$rMeans[2,]-mu[2]-bs[1,]-bs[3,]
    d1<-r*sum(dYb1^2)+parY$sv[1]
    d2<-r*sum(dYb2^2)+parY$sv[2]

    pA1=A[1]+n*r/2
    pA2=A[2]+n*r/2
    pB1=B[1]+d1/2
    pB2=B[2]+d2/2
    pA3=A[3]+n/2
    pB3=B[3]+sum((bs[1,])^2)/2

    c1<-pB1/rgamma(1,shape=pA1)    
    c2<-pB2/rgamma(1,shape=pA2)    
    c3<-pB3/rgamma(1,shape=pA3)
    if( identical (r, 1)) { c4 <- 0 } 
    else{
        pA4 <- A[4]+n
        pB4 <- B[4]+(sum(bs[2,]^2)+sum(bs[3,]^2))/2
        c4 <- 1/rgamma(1,shape=pA4,rate=pB4)
    }
    return(c(c1,c2,c3,c4))
}

################################

gamma.p<-function(theta,parY)
{
# computes full conditional of gamma(=mu,b) given theta and csY
 
    n<-parY$dim[1]
    r<-parY$dim[2]
  
    theta1r<-theta[1]/r;theta2r<-theta[2]/r

    l1<-sqrt(theta1r+theta[3]+theta[4])
    l3<-theta[3]/l1
    l2<-sqrt(theta2r+theta[3]+theta[4]-l3^2)
    N0<-rnorm(2)/sqrt(n)
    mu<-c(l1*N0[1],l3*N0[1]+l2*N0[2])+parY$mean
 
    lN1<-(N0[1]-N0[2]*l3/l2)/l1
    lN2<-N0[2]/l2

    t1<-theta1r+theta[4];t2<-theta2r+theta[4]
    g<-theta[3]/(t1*t2+theta[3]*(t1+t2))
 
    b<-sqrt(g*t1*t2)*rnorm(n)
    s1<-sqrt(theta1r*(theta[4]/t1))*rnorm(n)-(theta[4]/t1)*b
    s2<-sqrt(theta2r*(theta[4]/t2))*rnorm(n)-(theta[4]/t2)*b
 
    d1<-(t2/theta[3])*parY$dY[1,]
    d2<-(t1/theta[3])*parY$dY[2,]
    d12<-parY$dY[3,]
 
    b<-b+(g*theta[3])*(d1+d2)-theta[3]*(lN1+lN2)
    s1<-s1+(g*theta[4])*(d1+d12)-theta[4]*lN1
    s2<-s2+(g*theta[4])*(d2-d12)-theta[4]*lN2

    return(list(mu=mu,bs=rbind(b,s1,s2)))
 }

############################################################################