## com34a of september 2021,
#  re binormal data for a binormal prior 
#  parameters: (a,b); data: (x,y) 
#  what is Pr(a > b) for prior? 
#  what is Pr(a > b | data)? 

#  the prior: 
theta0 <- c(3.0,2.0)
dim(theta0) <- c(2,1) 
rho0 <- 0.88
Sigma0 <- rbind(c(1^2, 1*0.5*rho0), c(1*0.5*rho0,0.5^2))
print(Sigma0)

#  what is Pr(a > b) for prior? 
delta0 <- 3.0 - 2.0
tau0 <- sqrt(1^2 + 0.5^2 - 2*1*0.5*rho0)
#  delta0 + tau0 N > 0 : 
prob0 <- 1 - pnorm(-delta0/tau0) # 0.9035

#  the data observation process variance matrix: 
rho <- 0.33
Sigma <- rbind(c(1^2, 1*1*rho), c(1*1*rho,1^2))
print(Sigma) 

#  the data: 
xy <- c(2.77,2.33)
dim(xy) <- c(2,1) 
x <- xy[1]
y <- xy[2]

#  joint distribution for (a,b,x,y) is 4-normal, with variance matrix 
#  Sigma0, Sigma0
#  Sigma0, Sigma0 + Sigma 

# posterior is binormal: 
postmean <- theta0 + Sigma0 %*% solve(Sigma0 + Sigma) %*% (xy - theta0)
#  2.885
#  2.017 

postvariance <- Sigma0 - Sigma0 %*% solve(Sigma0+Sigma) %*% Sigma0

deltapost <- postmean[1] - postmean[2]
taupost <- sqrt(postvariance[1,1] + postvariance[2,2] - 2*postvariance[1,2])

#  a given data: N(2.8850, 0.7071^2)
#  b given data: N(2.0175, 0.4198^2) 
#  a-b given data: N(0.8675, 0.58884^2) 
prob <- 1 - pnorm(-deltapost/taupost) # 0.9298 

## lazy Bayesian, does (a,x) alone and (b,y) alone,
#  not caring about correlations in prior and in data: 

lazy1mean <- theta0[1] + (Sigma0[1,1]/(Sigma0[1,1]+Sigma[1,1])) * (x - theta0[1])
lazy1stdev <- sqrt(Sigma0[1,1] - Sigma0[1,1]^2/(Sigma0[1,1] + Sigma[1,1]))
# is N(2.8850, 0.7071^2) 

lazy2mean <- theta0[2] + (Sigma0[2,2]/(Sigma0[2,2]+Sigma[2,2])) * (y - theta0[2])
lazy2stdev <- sqrt(Sigma0[2,2] - Sigma0[2,2]^2/(Sigma0[2,2] + Sigma[2,2]))
# is N(2.066, 0.4472^2) 

print(postmean)
print(sqrt(postvariance))

print(c(lazy1mean,lazy1stdev))
print(c(lazy2mean,lazy2stdev))