#
# Ex.1a, Gordon et al., Prob.Prog., 2014 --- E[C1,C2] = <1/2, 1/2> ---
# 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.1a, Gordon et al., Prob.Prog., 2014 --- E[C1,C2] = <1/2, 1/2> ---")
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")
Trials <- matrix(c(C1, C2), byrow=TRUE, nrow=2) # Bernoulli trials
return(Trials)}
display <- function(Trials, title)
{nr <- length(Trials[1,])
cat("\neffective sample size nr = ", nr, "\n\n")
cat("E(C1) = ",mean(Trials[1,]), "\n")
cat("E(C2) = ",mean(Trials[2,]), "\n")
P <- table(Trials[1,],Trials[2,], dnn=c("C1", "C2"))/nr
cat("\nJPD P(C1, C2)\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("\nexpectation E(< C1, C2 > | C1=1 v C2=1) = < %f, %f > ",
EP[1], EP[2]))}
display(take_samples(n=150),
"JPD P(C1, C2)" )
##
## Ex.1a, Gordon et al., Prob.Prog., 2014 --- E[C1,C2] = <1/2, 1/2> ---
## R-Code by PCM 2015/04/20
##
## n = 150
## p = 0.5
##
## unfiltered trials
## C1 = 1 1 1 1 1 0 0 0 1 0 1 0 1 1 1 1 0 1 0 1 0 0 0 1 1 0 0 0 1 1 1 0 1 1 1 1 0 0 1 0 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 1 0 1 0 0 1 1 0 1 0 1 0 1 0 0 0 1 1 0 0 0 1 0 1 0 1 1 1 1 1 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 1 1 1 1 0 0 1 1 0 0 0 1 1 0 0
## C2 = 0 0 1 0 1 0 0 1 1 0 1 1 1 0 0 1 1 0 0 0 0 0 1 1 0 0 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 0 1 0 1 0 1 1 1 0 1 0 1 1 0 0 1 0 0 0 0 1 1 1 0 1 0 0 1 0 1 0 1 1 1 0 0 1 0 0 0 0 1 1 1 1 1 0 1 0 0 0 1 1 1 0 0 0 1 1 0 1 0 0 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 0 0 0 1 1 0 1 1 1 0 1 0 0 1 1 1 0 0 1 0
##
## effective sample size nr = 150
##
## E(C1) = 0.5
## E(C2) = 0.4733333
##
## JPD P(C1, C2)
## C2
## C1 0 1
## 0 0.2666667 0.2333333
## 1 0.2600000 0.2400000
##
## expectation E(< C1, C2 > | C1=1 v C2=1) = < 0.500000, 0.473333 >