#===============================================================================================# # The marginal posterior probability of (response B > response A) by normal approximation #===============================================================================================# # # MPP = function(v,alpha,obs) # (1) Choose a "v" value, a "alpha" vector and a "obs" vector as the number of response, # prior parameter vector and the observation, respectively. # (2) Note that both lengths of alpha observation is 2^v and here we assume that v is # lager than 2. # (3) For example, # # response <- 3 # prior.parameter <- c(5,98,63,7,42,7,7,7) # observation <- c(5,5,3,8,3,10,5,3) # MPP(response,prior.parameter,observation) # Output: # response A response B Pr(response B > response A) # [1,] 1 2 0.93692340 # [2,] 1 3 0.99999778 # [3,] 2 1 0.06307660 # [4,] 2 3 0.99812609 # [5,] 3 1 0.00000222 # [6,] 3 2 0.00187391 #==============================================================================================# MPP = function(v,alpha,obs) { A <- matrix(0,2^v,v) # create the parameter space (total of possible answers) for(j in 1:v){ A[,j] <- rep(c(rep(0,2^(v-j)),rep(1,2^(v-j))),2^(j-1)) } n <-(v-1)*v # create the numbers of all possible marginal posterior probability a1 <- numeric(n) a2 <- numeric(n) for(i in 1:v){ a1[((i-1)*(v-1)+1):((i-1)*(v-1)+(v-1))] <- rep(i,(v-1)) } a2[1:(v-1)] <- c(2:v) a2[((v-1)*(v-1)+1):((v-1)*(v-1)+(v-1))] <- c(1:(v-1)) for(i in 2:(v-1)){ a2[((v-1)*(i-1)+1):((v-1)*(i-1)+(v-1))] <- c(1:(i-1),(i+1):v) } AA <- cbind(A,alpha,obs) # combine the matrix A, the alpha and the observation V <- matrix(0,1,n) # create space that save numbers of all possible marginal posterior probability for(i in 1:n){ alpha0 <- sum(AA[,(v+1)]+AA[,(v+2)]) # compute A alphai <- AA[((AA[,a1[i]]!= AA[,a2[i]])& AA[,a1[i]]==1),][,(v+1)]+AA[((AA[,a1[i]]!=AA[,a2[i]])& AA[,a1[i]]==1),][,(v+2)] alphaj <- AA[((AA[,a1[i]]!= AA[,a2[i]])& AA[,a2[i]]==1),][,(v+1)]+AA[((AA[,a1[i]]!=AA[,a2[i]])& AA[,a2[i]]==1),][,(v+2)] E <- sum(alphai-alphaj)/alpha0 # compute B Var1 <- sum(alphai*alpha0-alphai^2) Var2 <- sum(alphaj*alpha0-alphaj^2) Var3 <- 2* sum(alphai%*%t(alphaj)) Var <- (Var1+Var2+Var3)/(alpha0^3+alpha0) # compute C V[1,i] <- pnorm(-E/sqrt(Var),0,1) # compute the approximation } tv <- round(V[1,],8) result <- cbind(a1,a2,tv) colnames(result) <- c("response A","response B","Pr(response B > response A)") return(result) }