#
# Ex.1b, Gordon et al., Prob.Prog., 2014 --- E[C1,C2 | c1 v c2] = <2/3, 2/3> ---
# R-Code by PCM 2015/04/20
#
library(latticeExtra) # needed for 3D-bivariate histogram (barchart)
## Loading required package: RColorBrewer
## Loading required package: lattice
library(knitr)
flips <- function(n, p) rbinom(n, 1, p) # generates n Bernoulli trials
or <- function(x, y) if(x | y) TRUE else FALSE # logical or function
expectation <- function(P) # expectation vector E(C1, C2)
{vecP <- as.vector(P)
M <- matrix(c(c(0,0),c(1,0),c(0,1),c(1,1)), nrow=2)
EP <- M %*% vecP} # matrix product M*P
take_samples <- function(n=20, p=0.5)
{cat("\nEx.1b, Gordon et al., Prob.Prog., 2014 --- E[C1,C2 | c1 v c2] = <2/3, 2/3> ---")
cat("\n R-Code by PCM 2015/04/20 \n")
cat("\nn = ", n)
cat("\np = ", p)
cat("\n\nunfiltered trials")
C1 <- flips(n, p) # n flips of coin1
cat("\nC1 = ", C1, "\n")
C2 <- flips(n, p) # n flips of coin2
cat("C2 = ", C2, "\n")
C1vC2 <- mapply(or, C1, C2) # C1=1 v C2=1
cat("\nC1 v C2 = \n", C1vC2, "\n")
FilteredTrials <- matrix(c(C1[C1vC2], C2[C1vC2]), byrow=TRUE, nrow=2) # conditioned Bernoulli trials
cat("\nfiltered trials\n")
print(FilteredTrials)
return(FilteredTrials)}
display <- function(Trials, title)
{nr <- length(Trials[1,])
cat("\neffective sample size nr = ", nr, "\n\n")
cat("E(C1 | C1 v C2) = ",mean(Trials[1,]), "\n")
cat("E(C2 | C1 v C2) = ",mean(Trials[2,]), "\n")
P <- table(Trials[1,],Trials[2,], dnn=c("C1", "C2"))/nr
cat("\nCPD P(C1, C2 | C1=1 v C2=2)\n")
print(P)
print(cloud(P,, panel.3d.cloud=panel.3dbars, col.facet='skyblue2',
main=title, xlab="C1", ylab="C2", zlab="P", zlim=c(0, 1),
xbase=0.5, ybase=0.5, scales=list(arrows=FALSE, col=1),
par.settings = list(axis.line = list(col = "transparent"))))
EP <- expectation(P)
cat(sprintf("\nconditional expectation E(< C1, C2 > | C1=1 v C2=1) = < %f, %f > ",
EP[1], EP[2]))}
display(take_samples(n=150),
"CPD P(C1, C2 | C1=1 v C2=1)" )
##
## Ex.1b, Gordon et al., Prob.Prog., 2014 --- E[C1,C2 | c1 v c2] = <2/3, 2/3> ---
## R-Code by PCM 2015/04/20
##
## n = 150
## p = 0.5
##
## unfiltered trials
## C1 = 0 0 1 1 1 1 0 1 0 0 1 1 0 0 0 1 1 1 0 0 1 0 1 0 0 0 1 0 1 1 1 0 0 1 0 1 0 0 1 0 1 0 1 1 0 0 1 1 0 1 1 1 1 1 1 0 0 1 1 0 1 1 1 1 0 1 0 1 0 1 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 0 0 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 0 1 0 0 1 0 1 1 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0
## C2 = 0 0 1 1 0 1 0 0 1 1 0 0 0 1 1 1 1 1 0 0 0 0 1 1 0 0 1 0 1 1 0 1 1 0 1 1 0 0 0 0 1 1 0 1 1 0 1 0 0 1 0 1 0 0 1 1 0 1 0 1 0 1 0 1 1 1 0 1 1 0 1 0 1 0 0 0 1 0 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 0 0 1 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 1 0 1 1 0 0 0 1 0 1 0 1 0 0 0 1 0 1 1 0 1 0 1 1 1 0 0 0 0 1 1 1 1 1 1 1
##
## C1 v C2 =
## FALSE FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##
## filtered trials
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 1 1 1 1 0 0 1 1 0 0 1 1
## [2,] 1 1 0 1 0 1 1 0 0 1 1 1 1
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 1 1 1 0 1 1 1 1 0 0 1
## [2,] 1 0 1 1 1 1 1 0 1 1 0
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 0 1 1 1 0 1 1 0 1 1 1
## [2,] 1 1 0 1 1 0 1 1 1 0 1
## [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,] 1 1 1 1 1 0 1 1 0 1 1
## [2,] 0 1 0 0 1 1 1 0 1 0 1
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57]
## [1,] 1 1 0 1 1 0 1 0 0 0 1
## [2,] 0 1 1 1 1 1 0 1 1 1 0
## [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,] 0 1 1 0 0 1 0 0 0 0 1
## [2,] 1 1 1 1 1 1 1 1 1 1 1
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79]
## [1,] 1 1 1 0 1 1 0 1 0 1 0
## [2,] 1 0 0 1 1 0 1 0 1 1 1
## [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] 1 0 0 1 0 1 0 1 1 0 1
## [2,] 0 1 1 0 1 1 1 1 0 1 1
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101]
## [1,] 1 0 0 1 0 1 0 1 1 0 0
## [2,] 0 1 1 0 1 0 1 1 1 1 1
## [,102] [,103] [,104] [,105] [,106] [,107] [,108] [,109] [,110] [,111]
## [1,] 1 1 1 0 0 1 0 0 1 0
## [2,] 1 1 0 1 1 1 1 1 1 1
##
## effective sample size nr = 111
##
## E(C1 | C1 v C2) = 0.6036036
## E(C2 | C1 v C2) = 0.7387387
##
## CPD P(C1, C2 | C1=1 v C2=2)
## C2
## C1 0 1
## 0 0.0000000 0.3963964
## 1 0.2612613 0.3423423
##
## conditional expectation E(< C1, C2 > | C1=1 v C2=1) = < 0.603604, 0.738739 >