In giuliomorina/DiceEnterprise: Bernoulli Factory and Dice Enterprise for Rational Functions
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width=7, fig.height=4,
cache=TRUE
)
require("DiceEnterprise")
require("pander")
require("ggplot2")
require("xtable")
set.seed(10,"L'Ecuyer-CMRG")
Example 1 - Toy Bernoulli Factory
$$
f(p) = \frac{\sqrt{2}p^3}{(\sqrt{2}-5)p^3+11p^2-9p+3}
$$
sample_size <- 1000 #Number of tosses
probs_vec <- c(0.01,0.1,0.25,0.5,0.75,0.9,0.99) #Vector of true probabilities of the p-coin
num_cores <- 2 #Number of cores used (parallel computing is not supported on Windows)
conf_level <- 0.95 #Confidence interval level
toss.coin <- function(n, p) {
sample(1:2, size = n, replace = TRUE, prob = c(p,1-p))
}
bf <- BernoulliFactory$new(f_1 = list(coeff = c(sqrt(2)), power = c(3)),
f_2 = list(coeff = c(-5,11,-9,3), power = c(3,2,1,0)))
bf.res <- function(toss.coin, bf) {
bf_res <- matrix(NA, nrow = 6, ncol = length(probs_vec))
rownames(bf_res) <- c("f(p)", "LowerCI", "Emp.f(p)", "UppCI", "Exp.TossesNoDoubling", "Exp.TossesDoubling")
colnames(bf_res) <- as.character(probs_vec)
for(i in 1:length(probs_vec)) {
sample_res <- bf$sample(n = sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE)
sample_res_double_time <- bf$sample(n = sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = TRUE)
#Get lower and upper CI
if(isTRUE(all.equal(sample_res[[1]],rep(1,sample_size)))) {
#All equal to 1 (Heads)
lower_ci <- 1
upper_ci <- 1
} else if(isTRUE(all.equal(sample_res[[1]],rep(2,sample_size)))) {
#All equal to 2 (Tails)
lower_ci <- 0
upper_ci <- 0
} else {
ci <- MultinomCI(table(sample_res[[1]]),conf.level = conf_level)
lower_ci <- ci[1,2]
upper_ci <- ci[1,3]
}
bf_res[1,i] <- round(bf$evaluate(p = probs_vec[i])[1],2)
bf_res[2,i] <- round(lower_ci,2)
bf_res[3,i] <- round(1-mean(sample_res[[1]]-1),2)
bf_res[4,i] <- round(upper_ci,2)
bf_res[5,i] <- round(mean(sample_res[[2]]),2)
bf_res[6,i] <- round(mean(sample_res_double_time[[2]]),2)
}
return(bf_res)
}
bf_res <- bf.res(toss.coin,bf)
pandoc.table(bf_res)
xtable(bf_res, align = c("l|",rep("c",ncol(bf_res))))
Example 2 - Efficiency of Monotonic CFTP
We consider the ladder
$$
\pi(p) \propto ((1-p)^4,1000p(1-p)^3,p^2(1-p)^2,500p^3(1-p),p^4)
$$
and show that increasing the degree of the ladder leads to better performance in terms of the expected number of tosses of the $p$-coin.
R <- c(1,1000,1,500,1)
max_increase_degree <- 2 #How much the degree of the original ladder is increased by
original_ladder <- DiceEnterprise$new(list(
list(R[1], "04"),
list(R[2], "13"),
list(R[3], "22"),
list(R[4], "31"),
list(R[5], "40")
))
#Sample and get the empirical number of tosses required
bf.res.eff <- function(toss.coin, original_ladder,max_increase_degree) {
#Prepare increased degree ladders
increased_ladders <- vector("list", length= max_increase_degree)
for(i in 1:max_increase_degree) {
increased_ladders[[i]] <- original_ladder$increase.degree(d=i)
}
#Prepare output
bf_res <- matrix(NA, nrow = max_increase_degree+1, ncol = length(probs_vec)+1)
rownames(bf_res) <- paste0("+",0:max_increase_degree)
colnames(bf_res) <- c(as.character(probs_vec),"Eff.Condition") #The last column is true if the efficiency condition is met
tot_rows <- length(probs_vec)*(max_increase_degree+1)
df_res <- data.frame(p=rep(NA,tot_rows),
degree = rep(NA, tot_rows),
empirical_tosses = rep(NA, tot_rows))
counter <- 1
#Compute
for(i in 1:length(probs_vec)) {
#Check if the result already exists
bf_res[1,i] <- mean(original_ladder$sample(n=sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE)$exp_rolls)
bf_res[1,ncol(bf_res)] <- as.numeric(original_ladder$get.efficiency.condition())
df_res[counter,] <- c(probs_vec[i],0,bf_res[1,i])
counter <- counter + 1
for(deg in 1:max_increase_degree) {
bf_res[deg+1,i] <- mean(increased_ladders[[deg]]$sample(n=sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE)$exp_rolls)
bf_res[deg+1,ncol(bf_res)] <- as.numeric(increased_ladders[[deg]]$get.efficiency.condition())
df_res[counter,] <- c(probs_vec[i],deg,bf_res[deg+1,i])
counter <- counter + 1
}
}
return(list(bf_res,df_res))
}
if(!file.exists(paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds"))) {
aux <- bf.res.eff(toss.coin,original_ladder,max_increase_degree)
#Save the result
saveRDS(list(aux=aux),file=paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds"))
} else {
aux <- readRDS(paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds"))$aux
}
bf_res <- aux[[1]]
print(bf_res)
pandoc.table(bf_res)
Example 3 - Dice Enterprise
$$
\pi(p_1,p_2,p_3) \propto \left(\sqrt{2}p_1^3,p_1^2p_3,\frac{1}{4}p_1p_2^2,2p_1p_2p_3,\frac{1}{2}p_1p_3^2,\frac{3}{4}p_2^2p_3\right)
$$
sample_size <- 1000 #Number of rolls
probs_mat <- matrix(c(1/5,1/4,11/20,
1/10,4/10,5/10), ncol = 3, byrow = TRUE) #Matrix of true probabilities of the die
num_cores <- 4 #Number of cores used (parallel computing is not supported on Windows)
conf_level <- 0.95 #Confidence interval level
roll_die <- function(n, probs_die) {
sample(1:3, size = n, replace = TRUE, prob = probs_die)
}
de <- DiceEnterprise$new(G=list(
list(sqrt(2),"300"),
list(1,"201"),
list(1/4,"120"),
list(2,"111"),
list(1/2,"102"),
list(3/4,"021")
))
de_res <- matrix(NA, nrow = 6, ncol = nrow(probs_mat))
rownames(de_res) <- c("f(p)", "LowerCI", "Emp.f(p)", "UppCI", "Exp.RollsNoDoubling", "Exp.RollsDoubling")
colnames(de_res) <- apply(probs_mat, 1, function(x) { paste0("(",paste0(x, collapse = ", "),")")} )
for(i in 1:nrow(probs_mat)) {
de_sample_double_time <- de$sample(n=sample_size, roll.fun = roll_die, double_time = TRUE, num_cores = num_cores, verbose = TRUE, probs_die = probs_mat[i,])
de_sample <- de$sample(n=sample_size, roll.fun = roll_die, double_time = FALSE, num_cores = num_cores, verbose = TRUE, probs_die = probs_mat[i,])
if(length(unique(de_sample[[1]])) > 1) {
#There are more outcomes
ci <- MultinomCI(table(de_sample[[1]]),conf.level = conf_level)
ci_aux <- matrix(NA, nrow = 2, ncol = 6) #6 possible outcomes of the die
for(j in 1:nrow(ci)) {
ci_aux[1,as.integer(rownames(ci)[j])] <- ci[j,2]
ci_aux[2,as.integer(rownames(ci)[j])] <- ci[j,3]
}
de_res[2,i] <- paste0("(",paste0(round(ci_aux[1,],2), collapse = ", "),")") #Lower confidence interval
de_res[4,i] <- paste0("(",paste0(round(ci_aux[2,],2), collapse = ", "),")") #Uppter confidence interval
}
de_res[1,i] <- paste0("(",paste0(round(de$evaluate(p = probs_mat[i,]),2), collapse = ", "),")")
de_res[3,i] <- paste0("(",paste0(round(table(de_sample[[1]])/sample_size,2), collapse = ", "),")")
de_res[5,i] <- round(mean(de_sample[[2]]),2)
de_res[6,i] <- round(mean(de_sample_double_time[[2]]),2)
}
pandoc.table(de_res,split.tables=Inf)
Example 4 - Augmenting ladder (Bernoulli factory)
We consider a ladder of the following form, where the true (unknown) $p$ is $1/2$:
$$
\pi(p) \propto {p^{2a},p^a(1-p)^a,(1-p)^{2a}}
$$
A naive approach similar to the 2-coin algorithm would be with probability $1/3$, toss a $p$-coin $2a$ times and output $1$ if all heads; with prob. $1/3$ toss a $p$-coin $2a$ times and output $1$ if the first $a$ are heads and the last $a$ tosses are tails; with prob. $1/3$ toss a $p$-coin $2a$ times and output $1$ if all tails. This would require on average $\mathbb{E}_{p=1/2}[N] = 2^{2a}$.
rm(list=ls())
require("DiceEnterprise")
p <- 1/2 #True value of p
a <- 10 #Exponent
n <- 1000 #Size sample CFTP
num_cores <- 4 #Number of cores used
initial_seed <- 17 #Seed for RNG
increase_degree <- 260 #How much the degree is increased
de_list <- vector("list", length = increase_degree+1) #Store all the de objects
set.seed(initial_seed, "L'Ecuyer-CMRG")
#Construct the ladder
cat("Constructing the dice enterprise \n")
de <- DiceEnterprise$new(G=list(
list(1,matrix(c(2*a,0),nrow=1,ncol=2)),
list(1,matrix(c(a,a),nrow=1,ncol=2)),
list(1,matrix(c(0,2*a),nrow=1,ncol=2))
))
de_list[[1]] <- de
cat("Constructing the augmented dice enterprises \n")
for(i in 1:increase_degree) {
#Check if it already exists
if(file.exists(paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds"))) {
cat("Restoring previous work for degree +",i," \n")
prev_work <- readRDS(paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds"))
de_list[[i+1]] <- prev_work$de_aug
} else {
stop("WHAT?!")
cat("Constructing the augmented dice enterprises for +",i," \n")
if(i == 1) {
de_list[[i+1]] <- de$increase.degree(1)
} else {
de_list[[i+1]] <- de_list[[i-1]]$increase.degree(1)
}
#Save object
saveRDS(list(de_aug=de_list[[i+1]],i=i), file = paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds"))
}
}
#For each ladder get empirical number of tosses and bound
cat("Computing empirical number of tosses and bound \n")
emp_tosses <- matrix(NA, nrow = n, ncol = increase_degree+1)
emp_res <- matrix(NA, nrow = n, ncol = increase_degree+1)
colnames(emp_tosses) <- paste0("+",0:increase_degree)
colnames(emp_res) <- paste0("+",0:increase_degree)
for(i in 0:increase_degree) {
#Check if it already exists
if(file.exists(paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds"))) {
cat("Restoring previous work for degree +",i," \n")
prev_work <- readRDS(paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds"))
emp_tosses[,i+1] <- prev_work$emp_tosses
emp_res[,i+1] <- prev_work$emp_res
} else {
cat("Getting empirical tosses for +",i," \n")
aux <- de_list[[i+1]]$sample(n=n, true_p = c(p,1-p), num_cores = num_cores, verbose = TRUE)
emp_tosses[,i+1] <- aux[[2]]
emp_res[,i+1] <- aux[[1]]
#Save object
saveRDS(list(emp_tosses=aux[[2]],emp_res=aux[[1]]), file = paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds"))
}
}
#Prepare output
res_output <- matrix(NA, nrow = 2, ncol = increase_degree+1)
colnames(res_output) <- paste0("+",0:increase_degree)
rownames(res_output) <- c("Exp.Tosses","Efficiency.Condition")
res_output[1,] <- apply(emp_tosses,2,mean)
res_output[2,] <- sapply(de_list, function(x) {x$get.efficiency.condition()})
print(res_output)
Example 5 - Augmenting ladder (Dice enterprise)
We consider a ladder of the following form, where the true (unknown) $\boldsymbol{p}=(p_1,p_2,p_3)$ is $(1/3,1/3,1/3)$:
$$
\pi(\boldsymbol{p}) \propto {p_1^ap_2^ap_3^a,p_1^{3a},p_2^{3a},p_3^{3a}}
$$
rm(list=ls())
require("DiceEnterprise")
setwd("~/Dropbox/OxWaSP/Dice Enterprise/R/DiceEnterprise/vignettes") #Change accordingly
true_p <- c(1/3,1/3,1/3) #True value of p
a <- 5 #Exponent
n <- 1000 #Size sample CFTP
num_cores <- 4 #Number of cores used
initial_seed <- 7117 #Seed for RNG
increase_degree <- 78 #How much the degree is increased
de_list <- vector("list", length = increase_degree+1) #Store all the de objects
set.seed(initial_seed, "L'Ecuyer-CMRG")
#Construct the ladder
cat("Constructing the dice enterprise \n")
de <- DiceEnterprise$new(G=list(
list(1,matrix(c(3*a,0,0),nrow=1,ncol=3)),
list(1,matrix(c(a,a,a),nrow=1,ncol=3)),
list(1,matrix(c(0,0,3*a),nrow=1,ncol=3)),
list(1,matrix(c(0,3*a,0),nrow=1,ncol=3))
))
de_list[[1]] <- de
cat("Constructing the augmented dice enterprises \n")
for(i in 1:increase_degree) {
#Check if it already exists
if(file.exists(paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds"))) {
cat("Restoring previous work for degree +",i," \n")
prev_work <- readRDS(paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds"))
de_list[[i+1]] <- prev_work$de_aug
} else {
cat("Constructing the augmented dice enterprises for +",i," \n")
if(i == 1) {
de_list[[i+1]] <- de$increase.degree(1)
} else {
de_list[[i+1]] <- de_list[[i-1]]$increase.degree(1)
}
#Save object
saveRDS(list(de_aug=de_list[[i+1]],i=i), file = paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds"))
}
}
#For each ladder get empirical number of tosses
cat("Computing empirical number of tosses \n")
emp_tosses <- matrix(NA, nrow = n, ncol = increase_degree+1)
emp_res <- matrix(NA, nrow = n, ncol = increase_degree+1)
colnames(emp_tosses) <- paste0("+",0:increase_degree)
colnames(emp_res) <- paste0("+",0:increase_degree)
for(i in 0:increase_degree) {
#Check if it already exists
if(file.exists(paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds"))) {
cat("Restoring previous work for degree +",i," \n")
prev_work <- readRDS(paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds"))
emp_tosses[,i+1] <- prev_work$emp_tosses
emp_res[,i+1] <- prev_work$emp_res
} else {
cat("Getting empirical tosses for +",i," \n")
aux <- de_list[[i+1]]$sample(n=n, true_p = true_p, num_cores = num_cores, verbose = TRUE)
emp_tosses[,i+1] <- aux[[2]]
emp_res[,i+1] <- aux[[1]]
#Save object
saveRDS(list(emp_tosses=aux[[2]],emp_res=aux[[1]]), file = paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds"))
}
}
apply(emp_tosses,2,mean)
#Prepare output
res_output <- matrix(NA, nrow = 1, ncol = increase_degree+1)
colnames(res_output) <- paste0("+",0:increase_degree)
res_output[1,] <- apply(emp_tosses,2,mean)
print(res_output)
giuliomorina/DiceEnterprise documentation built on May 15, 2019, 12:59 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=7, fig.height=4, cache=TRUE ) require("DiceEnterprise") require("pander") require("ggplot2") require("xtable") set.seed(10,"L'Ecuyer-CMRG")
Example 1 - Toy Bernoulli Factory
$$ f(p) = \frac{\sqrt{2}p^3}{(\sqrt{2}-5)p^3+11p^2-9p+3} $$
sample_size <- 1000 #Number of tosses probs_vec <- c(0.01,0.1,0.25,0.5,0.75,0.9,0.99) #Vector of true probabilities of the p-coin num_cores <- 2 #Number of cores used (parallel computing is not supported on Windows) conf_level <- 0.95 #Confidence interval level
toss.coin <- function(n, p) { sample(1:2, size = n, replace = TRUE, prob = c(p,1-p)) } bf <- BernoulliFactory$new(f_1 = list(coeff = c(sqrt(2)), power = c(3)), f_2 = list(coeff = c(-5,11,-9,3), power = c(3,2,1,0))) bf.res <- function(toss.coin, bf) { bf_res <- matrix(NA, nrow = 6, ncol = length(probs_vec)) rownames(bf_res) <- c("f(p)", "LowerCI", "Emp.f(p)", "UppCI", "Exp.TossesNoDoubling", "Exp.TossesDoubling") colnames(bf_res) <- as.character(probs_vec) for(i in 1:length(probs_vec)) { sample_res <- bf$sample(n = sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE) sample_res_double_time <- bf$sample(n = sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = TRUE) #Get lower and upper CI if(isTRUE(all.equal(sample_res[[1]],rep(1,sample_size)))) { #All equal to 1 (Heads) lower_ci <- 1 upper_ci <- 1 } else if(isTRUE(all.equal(sample_res[[1]],rep(2,sample_size)))) { #All equal to 2 (Tails) lower_ci <- 0 upper_ci <- 0 } else { ci <- MultinomCI(table(sample_res[[1]]),conf.level = conf_level) lower_ci <- ci[1,2] upper_ci <- ci[1,3] } bf_res[1,i] <- round(bf$evaluate(p = probs_vec[i])[1],2) bf_res[2,i] <- round(lower_ci,2) bf_res[3,i] <- round(1-mean(sample_res[[1]]-1),2) bf_res[4,i] <- round(upper_ci,2) bf_res[5,i] <- round(mean(sample_res[[2]]),2) bf_res[6,i] <- round(mean(sample_res_double_time[[2]]),2) } return(bf_res) } bf_res <- bf.res(toss.coin,bf) pandoc.table(bf_res)
xtable(bf_res, align = c("l|",rep("c",ncol(bf_res))))
Example 2 - Efficiency of Monotonic CFTP
We consider the ladder $$ \pi(p) \propto ((1-p)^4,1000p(1-p)^3,p^2(1-p)^2,500p^3(1-p),p^4) $$ and show that increasing the degree of the ladder leads to better performance in terms of the expected number of tosses of the $p$-coin.
R <- c(1,1000,1,500,1) max_increase_degree <- 2 #How much the degree of the original ladder is increased by original_ladder <- DiceEnterprise$new(list( list(R[1], "04"), list(R[2], "13"), list(R[3], "22"), list(R[4], "31"), list(R[5], "40") )) #Sample and get the empirical number of tosses required bf.res.eff <- function(toss.coin, original_ladder,max_increase_degree) { #Prepare increased degree ladders increased_ladders <- vector("list", length= max_increase_degree) for(i in 1:max_increase_degree) { increased_ladders[[i]] <- original_ladder$increase.degree(d=i) } #Prepare output bf_res <- matrix(NA, nrow = max_increase_degree+1, ncol = length(probs_vec)+1) rownames(bf_res) <- paste0("+",0:max_increase_degree) colnames(bf_res) <- c(as.character(probs_vec),"Eff.Condition") #The last column is true if the efficiency condition is met tot_rows <- length(probs_vec)*(max_increase_degree+1) df_res <- data.frame(p=rep(NA,tot_rows), degree = rep(NA, tot_rows), empirical_tosses = rep(NA, tot_rows)) counter <- 1 #Compute for(i in 1:length(probs_vec)) { #Check if the result already exists bf_res[1,i] <- mean(original_ladder$sample(n=sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE)$exp_rolls) bf_res[1,ncol(bf_res)] <- as.numeric(original_ladder$get.efficiency.condition()) df_res[counter,] <- c(probs_vec[i],0,bf_res[1,i]) counter <- counter + 1 for(deg in 1:max_increase_degree) { bf_res[deg+1,i] <- mean(increased_ladders[[deg]]$sample(n=sample_size, roll.fun = toss.coin, p = probs_vec[i], num_cores = num_cores, verbose = TRUE, double_time = FALSE)$exp_rolls) bf_res[deg+1,ncol(bf_res)] <- as.numeric(increased_ladders[[deg]]$get.efficiency.condition()) df_res[counter,] <- c(probs_vec[i],deg,bf_res[deg+1,i]) counter <- counter + 1 } } return(list(bf_res,df_res)) } if(!file.exists(paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds"))) { aux <- bf.res.eff(toss.coin,original_ladder,max_increase_degree) #Save the result saveRDS(list(aux=aux),file=paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds")) } else { aux <- readRDS(paste0("saved_data/efficiency_mono_cftp_example_n_",sample_size,".rds"))$aux } bf_res <- aux[[1]] print(bf_res) pandoc.table(bf_res)
Example 3 - Dice Enterprise
$$ \pi(p_1,p_2,p_3) \propto \left(\sqrt{2}p_1^3,p_1^2p_3,\frac{1}{4}p_1p_2^2,2p_1p_2p_3,\frac{1}{2}p_1p_3^2,\frac{3}{4}p_2^2p_3\right) $$
sample_size <- 1000 #Number of rolls probs_mat <- matrix(c(1/5,1/4,11/20, 1/10,4/10,5/10), ncol = 3, byrow = TRUE) #Matrix of true probabilities of the die num_cores <- 4 #Number of cores used (parallel computing is not supported on Windows) conf_level <- 0.95 #Confidence interval level
roll_die <- function(n, probs_die) { sample(1:3, size = n, replace = TRUE, prob = probs_die) } de <- DiceEnterprise$new(G=list( list(sqrt(2),"300"), list(1,"201"), list(1/4,"120"), list(2,"111"), list(1/2,"102"), list(3/4,"021") )) de_res <- matrix(NA, nrow = 6, ncol = nrow(probs_mat)) rownames(de_res) <- c("f(p)", "LowerCI", "Emp.f(p)", "UppCI", "Exp.RollsNoDoubling", "Exp.RollsDoubling") colnames(de_res) <- apply(probs_mat, 1, function(x) { paste0("(",paste0(x, collapse = ", "),")")} ) for(i in 1:nrow(probs_mat)) { de_sample_double_time <- de$sample(n=sample_size, roll.fun = roll_die, double_time = TRUE, num_cores = num_cores, verbose = TRUE, probs_die = probs_mat[i,]) de_sample <- de$sample(n=sample_size, roll.fun = roll_die, double_time = FALSE, num_cores = num_cores, verbose = TRUE, probs_die = probs_mat[i,]) if(length(unique(de_sample[[1]])) > 1) { #There are more outcomes ci <- MultinomCI(table(de_sample[[1]]),conf.level = conf_level) ci_aux <- matrix(NA, nrow = 2, ncol = 6) #6 possible outcomes of the die for(j in 1:nrow(ci)) { ci_aux[1,as.integer(rownames(ci)[j])] <- ci[j,2] ci_aux[2,as.integer(rownames(ci)[j])] <- ci[j,3] } de_res[2,i] <- paste0("(",paste0(round(ci_aux[1,],2), collapse = ", "),")") #Lower confidence interval de_res[4,i] <- paste0("(",paste0(round(ci_aux[2,],2), collapse = ", "),")") #Uppter confidence interval } de_res[1,i] <- paste0("(",paste0(round(de$evaluate(p = probs_mat[i,]),2), collapse = ", "),")") de_res[3,i] <- paste0("(",paste0(round(table(de_sample[[1]])/sample_size,2), collapse = ", "),")") de_res[5,i] <- round(mean(de_sample[[2]]),2) de_res[6,i] <- round(mean(de_sample_double_time[[2]]),2) } pandoc.table(de_res,split.tables=Inf)
Example 4 - Augmenting ladder (Bernoulli factory)
We consider a ladder of the following form, where the true (unknown) $p$ is $1/2$: $$ \pi(p) \propto {p^{2a},p^a(1-p)^a,(1-p)^{2a}} $$ A naive approach similar to the 2-coin algorithm would be with probability $1/3$, toss a $p$-coin $2a$ times and output $1$ if all heads; with prob. $1/3$ toss a $p$-coin $2a$ times and output $1$ if the first $a$ are heads and the last $a$ tosses are tails; with prob. $1/3$ toss a $p$-coin $2a$ times and output $1$ if all tails. This would require on average $\mathbb{E}_{p=1/2}[N] = 2^{2a}$.
rm(list=ls()) require("DiceEnterprise") p <- 1/2 #True value of p a <- 10 #Exponent n <- 1000 #Size sample CFTP num_cores <- 4 #Number of cores used initial_seed <- 17 #Seed for RNG increase_degree <- 260 #How much the degree is increased de_list <- vector("list", length = increase_degree+1) #Store all the de objects set.seed(initial_seed, "L'Ecuyer-CMRG") #Construct the ladder cat("Constructing the dice enterprise \n") de <- DiceEnterprise$new(G=list( list(1,matrix(c(2*a,0),nrow=1,ncol=2)), list(1,matrix(c(a,a),nrow=1,ncol=2)), list(1,matrix(c(0,2*a),nrow=1,ncol=2)) )) de_list[[1]] <- de cat("Constructing the augmented dice enterprises \n") for(i in 1:increase_degree) { #Check if it already exists if(file.exists(paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds"))) { cat("Restoring previous work for degree +",i," \n") prev_work <- readRDS(paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds")) de_list[[i+1]] <- prev_work$de_aug } else { stop("WHAT?!") cat("Constructing the augmented dice enterprises for +",i," \n") if(i == 1) { de_list[[i+1]] <- de$increase.degree(1) } else { de_list[[i+1]] <- de_list[[i-1]]$increase.degree(1) } #Save object saveRDS(list(de_aug=de_list[[i+1]],i=i), file = paste0("efficiency_uni_data/ladder_a_",a,"_degree_",i,".rds")) } } #For each ladder get empirical number of tosses and bound cat("Computing empirical number of tosses and bound \n") emp_tosses <- matrix(NA, nrow = n, ncol = increase_degree+1) emp_res <- matrix(NA, nrow = n, ncol = increase_degree+1) colnames(emp_tosses) <- paste0("+",0:increase_degree) colnames(emp_res) <- paste0("+",0:increase_degree) for(i in 0:increase_degree) { #Check if it already exists if(file.exists(paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds"))) { cat("Restoring previous work for degree +",i," \n") prev_work <- readRDS(paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds")) emp_tosses[,i+1] <- prev_work$emp_tosses emp_res[,i+1] <- prev_work$emp_res } else { cat("Getting empirical tosses for +",i," \n") aux <- de_list[[i+1]]$sample(n=n, true_p = c(p,1-p), num_cores = num_cores, verbose = TRUE) emp_tosses[,i+1] <- aux[[2]] emp_res[,i+1] <- aux[[1]] #Save object saveRDS(list(emp_tosses=aux[[2]],emp_res=aux[[1]]), file = paste0("efficiency_uni_data/n_",n,"_a_",a,"_degree_",i,".rds")) } }
#Prepare output res_output <- matrix(NA, nrow = 2, ncol = increase_degree+1) colnames(res_output) <- paste0("+",0:increase_degree) rownames(res_output) <- c("Exp.Tosses","Efficiency.Condition") res_output[1,] <- apply(emp_tosses,2,mean) res_output[2,] <- sapply(de_list, function(x) {x$get.efficiency.condition()}) print(res_output)
Example 5 - Augmenting ladder (Dice enterprise)
We consider a ladder of the following form, where the true (unknown) $\boldsymbol{p}=(p_1,p_2,p_3)$ is $(1/3,1/3,1/3)$: $$ \pi(\boldsymbol{p}) \propto {p_1^ap_2^ap_3^a,p_1^{3a},p_2^{3a},p_3^{3a}} $$
rm(list=ls()) require("DiceEnterprise") setwd("~/Dropbox/OxWaSP/Dice Enterprise/R/DiceEnterprise/vignettes") #Change accordingly true_p <- c(1/3,1/3,1/3) #True value of p a <- 5 #Exponent n <- 1000 #Size sample CFTP num_cores <- 4 #Number of cores used initial_seed <- 7117 #Seed for RNG increase_degree <- 78 #How much the degree is increased de_list <- vector("list", length = increase_degree+1) #Store all the de objects set.seed(initial_seed, "L'Ecuyer-CMRG") #Construct the ladder cat("Constructing the dice enterprise \n") de <- DiceEnterprise$new(G=list( list(1,matrix(c(3*a,0,0),nrow=1,ncol=3)), list(1,matrix(c(a,a,a),nrow=1,ncol=3)), list(1,matrix(c(0,0,3*a),nrow=1,ncol=3)), list(1,matrix(c(0,3*a,0),nrow=1,ncol=3)) )) de_list[[1]] <- de cat("Constructing the augmented dice enterprises \n") for(i in 1:increase_degree) { #Check if it already exists if(file.exists(paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds"))) { cat("Restoring previous work for degree +",i," \n") prev_work <- readRDS(paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds")) de_list[[i+1]] <- prev_work$de_aug } else { cat("Constructing the augmented dice enterprises for +",i," \n") if(i == 1) { de_list[[i+1]] <- de$increase.degree(1) } else { de_list[[i+1]] <- de_list[[i-1]]$increase.degree(1) } #Save object saveRDS(list(de_aug=de_list[[i+1]],i=i), file = paste0("efficiency_multi_data/ladder_n_",n,"_a_",a,"_degree_",i,".rds")) } } #For each ladder get empirical number of tosses cat("Computing empirical number of tosses \n") emp_tosses <- matrix(NA, nrow = n, ncol = increase_degree+1) emp_res <- matrix(NA, nrow = n, ncol = increase_degree+1) colnames(emp_tosses) <- paste0("+",0:increase_degree) colnames(emp_res) <- paste0("+",0:increase_degree) for(i in 0:increase_degree) { #Check if it already exists if(file.exists(paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds"))) { cat("Restoring previous work for degree +",i," \n") prev_work <- readRDS(paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds")) emp_tosses[,i+1] <- prev_work$emp_tosses emp_res[,i+1] <- prev_work$emp_res } else { cat("Getting empirical tosses for +",i," \n") aux <- de_list[[i+1]]$sample(n=n, true_p = true_p, num_cores = num_cores, verbose = TRUE) emp_tosses[,i+1] <- aux[[2]] emp_res[,i+1] <- aux[[1]] #Save object saveRDS(list(emp_tosses=aux[[2]],emp_res=aux[[1]]), file = paste0("efficiency_multi_data/n_",n,"_a_",a,"_degree_",i,".rds")) } } apply(emp_tosses,2,mean)
#Prepare output res_output <- matrix(NA, nrow = 1, ncol = increase_degree+1) colnames(res_output) <- paste0("+",0:increase_degree) res_output[1,] <- apply(emp_tosses,2,mean) print(res_output)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.