Technical_Appendix.Rmd

---
title: "Technical Appendix"
output: html_document
date: "2024-07-21"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 

# Technical Appendix

## R code for Chapter 1

There is no code for chapter 1.

## R code for Chapter 2



```{r eval=F}

# Load the required packages
library(bayesAB)
library(dplyr)
library(ggplot2)

# Set seed for reproducibility
set.seed(123)

# Simulate data for group A
n_A <- 1000
conversion_rate_A <- 0.40
data_A <- rbinom(n_A, 1, conversion_rate_A)

# Define the range of conversion rates for group B
conversion_rate_B_values <- seq(0, 1, by = 0.01)
n_B <- 1000

# Initialize vectors to store results
p_values <- A_est <- B_est<-  
  numeric(length(conversion_rate_B_values))
bayesian_probabilities <- 
  numeric(length(conversion_rate_B_values))

# Informative priors
prior_alpha_A <- 10  # Prior successes 
prior_beta_A <- 90   # Prior failures 


# Loop over each conversion rate for group B
for (i in seq_along(conversion_rate_B_values)) {
  conversion_rate_B <- conversion_rate_B_values[i]
  
  # Generate data for group B
  data_B <- rbinom(n_B, 1, conversion_rate_B)
  
  # Frequentist A/B test
  freq_test <- prop.test(c(sum(data_A),
                           sum(data_B)), 
                         c(n_A, n_B))
  p_values[i] <- freq_test$p.value
  A_est[i] <- freq_test$estimate[[1]]
  B_est[i] <- freq_test$estimate[[2]]
  
  # Bayesian A/B test
  bayes_test <- bayesAB::bayesTest(
    data_A,
    data_B,
    priors = c('alpha' = prior_alpha_A,
               'beta' = prior_beta_A),
    distribution = 'bernoulli'
  )
  
  # Extract the probability that B is better than A
  bayesian_probabilities[i] <- 
    max(bayes_test$posteriors$Probability$B)
}

# Create a data frame for plotting
results_df <- data.frame(
  conversion_rate_B = 
    conversion_rate_B_values,
  p_value = p_values,
  bayesian_probability =
    bayesian_probabilities,
  A_values = A_est, 
  B_values = B_est
)


ggplot(results_df) +
  geom_line(aes(x = conversion_rate_B, 
                y = p_value), color = 'darkgrey') +
  geom_line(aes(x = conversion_rate_B, 
                y = bayesian_probability), 
            color = 'black') +
  geom_hline(yintercept = 0.4, 
             linetype = 'dashed', 
             color = 'lightgrey') +
  geom_vline(xintercept = 0.4, 
             linetype = 'dashed', 
             color = 'lightgrey') +
  labs(title = " ", 
#'Bayesian A/B Probabilities (black) versus 
#'Frequentist p-values (grey)',
       x = 'Conversion Rate B',
       y = 'Probability B > A & p-value') +
       theme_bw() 


```



```{r eval=F}

library(ggplot2)
library(tidyr)
library(reshape2)
library(tidyverse)
library(kableExtra)

stats <- data.frame()
hypt <-  vector()
for (j in 1:2000) {
  hypt[j] <- i <-  0
  for (i in 1:100) {
    k <- j * 5
    A <- rnorm(k, 110000, 300000)
    B <- rnorm(k, 100000, 300000)
    if (t.test(x = A, y = B, 
               paired = TRUE)$statistic < 0)
      hypt[j] <- hypt[j] + 1
    stats[i, j] <- t.test(x = A, y = B,
                          paired = TRUE)$p.value
  }
}

meanp <- apply(stats, 2, mean)
mp <- data.frame(meanp, seq(5, 10000, 5))
colnames(mp) <- c("p", "sample_size")

hyp <- cbind(hypt, seq(5, 10000, 5)) %>% 
  as.data.frame() 
colnames(hyp) <- c("pct_false", "sample_size")
to_graph <- 
  inner_join(mp, hyp,
             by="sample_size") %>%  
  mutate(pct_false = pct_false/100)

long_graph <- pivot_longer(to_graph, 
                           cols = c("p","pct_false"))
long_graph <- long_graph %>% 
  mutate(name = ifelse(name == "p", 
                       "p-value", 
                       "Wrong Decisions"))
ggplot(data=long_graph, aes(x=sample_size, 
                            y=value, 
                            colour=name)) +
  geom_smooth() +
    scale_colour_grey(start = 0, end = .6) + 
  theme_bw() +
#  ggtitle("More Information Lowers the Proportion of Wrong Decisions") +
  xlab("Sample Size") + 
  ylab("Proportion") + 
  theme(legend.title=element_blank())


```


## R code for Chapter 3



```{r eval=F}
library(bayesAB)
library(tidyverse)

# Set the lambda values for datasets A and B
lambda_A <- 45
lambda_B <- 46

# Define the size of each dataset
size_A <- 1000  
# Number of data points in dataset A
size_B <- 1000 
# Number of data points in dataset B

# Generate Dataset A with Poisson distribution
dataset_A <- rpois(size_A, lambda_A)

# Generate Dataset B with Poisson distribution
dataset_B <- rpois(size_B, lambda_B)

# Optionally shuffle data to randomize the order
set.seed(42)  # Setting seed for reproducibility
dataset_A <- sample(dataset_A)
dataset_B <- sample(dataset_B)


# Scenario 1: Weak prior 
test_weak_prior <- bayesTest(
  A = dataset_A,
  B = dataset_B,
  priors = c("shape" = 1, "rate" = 1),
  dist = 'poisson'
)

plot(test_weak_prior)
```

## R code for Chapter 4


```{r eval=F}
library(tidyverse)
library(ggplot2)

# Define the parameters for the two normal distributions
mean1 <- 600
sd1 <- 200
mean2 <- 400
sd2 <- 200


q_B <- qnorm(seq(0,1,.01), 
             mean =  mean2, sd = sd2)
sq <- seq(0,1,.01)
q_B <- cbind(sq,q_B) %>% 
  as.data.frame()  

# prob=.85 implies .15 in tail  



q_A <- qnorm(seq(0,1,.01), 
             mean =  mean1, sd = sd1)
q_A <- cbind(sq,q_A) %>% 
  as.data.frame()

## prob=.98 implies .02 in tail


# Create a sequence of x values
x <- seq(00, 1500, length.out = 1000)

# Calculate the density values for both distributions
density1 <- dnorm(x, 
                  mean = mean1,
                  sd = sd1)
density2 <- dnorm(x, 
                  mean = mean2, 
                  sd = sd2)

# Create a data frame to hold the values
df <- data.frame(x = c(x, x),
                 density = 
                   c(density1, density2),
                 group = 
                   factor(rep(c("Choice A", 
                                "Choice B"),
                                    each = length(x))))

# Plot the density functions
ggplot(df, aes(x = x, y = density, color = group)) +
  geom_line(size = 1) +
  scale_color_manual(values = c("grey", "black")) +
  theme_minimal() +  
  geom_vline(xintercept = 800,
             linetype = "dashed") +
  labs(title = 
         "Density Functions of Two Posterior Distributions",
       x = "Unit sales",
       y = "Density") +
  theme(legend.title = element_blank()) + xlim(0,1500)  

```


```{r eval=F}

library(tidyverse)
library(ggplot2)

# Define the profit function
profit_function <- function(units) {
  profit <- 4 * units - 1000
  return(profit)
}

# Define the density functions for the two distributions
density_1 <- function(x) {
  dnorm(x, mean = 600, sd = 200)
}

density_2 <- function(x) {
  dnorm(x, mean = 400, sd = 200)
}

# Generate data for plotting
units_sold <- seq(0, 1000, by = 1)
profit_1 <- 
  density_1(units_sold) * 
  profit_function(units_sold) * units_sold
profit_2 <- 
  density_2(units_sold) * 
  profit_function(units_sold) * units_sold

# Create a data frame for ggplot
data <- data.frame(
  Units = rep(units_sold, 2),
  Profit = c(profit_1, profit_2),
  Choice = rep(c("Choice A", "Choice B"), 
               each = length(units_sold))
)

# Plot the data
ggplot(data, aes(x = Units, y = Profit, color = Choice)) +
  geom_line(lwd=1) +
  scale_color_manual(values = c("grey50", "grey80")) +
  labs(
    title = "Expected Total Profit for Number of Units Sold",
    x = "Number of Units Sold",
    y = "Expected Total Profit",
    caption = "Profit maximizing choice"
  ) +
  theme_minimal()



```


## R code for Chapter 5

```{r eval=F}

library(tidyverse)
library(bayesAB)

set.seed(123456)

A_MW_IMDB <- rpois(10000, 2)
B_MotB_IMDB <- rpois(10000, 3)

AB_IMDB <- bayesTest(A_MW_IMDB, B_MotB_IMDB, priors =
                       c('shape' = 4, 'rate' = 4),
                     distribution = 'poisson')

A_MW_RT <- rpois(10000, 3)
B_MotB_RT <- rpois(10000, 2)

AB_RT <- bayesTest(
  A_MW_RT, 
  B_MotB_RT, 
  priors = c('shape' = 4, 'rate' = 4), 
  distribution = 'poisson')



AB <- combine(AB_IMDB, AB_RT, f = `+`)

plot(AB)

```



## R code for Chapter 6



```{r eval=F}

library(bayesAB)
library(tidyverse)
library(readr)

coin <- 
  read_csv("DATABASE_xisBTC_yisETH_coin.csv") %>% 
  select(btc=High.x, eth=High.y, date=Date)

# Example price series for Asset A and Asset B
price_btc <-  coin$btc
#  c(100, 102, 105, 110, 108, 115) # Example prices for Asset A
price_eth <- coin$eth
# c(50, 52, 55, 60, 58, 65)      # Example prices for Asset B

compute_incremental_roi <- function(prices) {
  # Calculate ROI between each successive data point
  roi <- diff(prices) / head(prices, -1)
  return(roi)
}

# Compute ROI for both series
roi_a <- compute_incremental_roi(price_btc)
roi_b <- compute_incremental_roi(price_eth)


df_btc <- data.frame(roi_a,"BTC")
colnames(df_btc) <- c("ROI", "Cryptocurrency")
df_eth <- data.frame(roi_b,"ETH")
colnames(df_eth) <- c("ROI", "Cryptocurrency")
df_coin <- rbind(df_btc,df_eth)

df_coin %>% ggplot(aes(ROI,color=Cryptocurrency))+
  geom_density(lwd=1)+theme_minimal()

# Perform the Bayes A/B test
ab_test_result <- bayesTest(
  A = roi_a,
  B = roi_b,
  priors = c("mu" = 0, 
             "lambda" = 1, 
             "alpha" = 2, 
             "beta" = .05),
  n_samples = 1e+05,
  dist = 'normal'
)

plot(ab_test_result)
#summary(ab_test_result)

```



```{r eval=F}

library(tidyverse)
library(reshape2)
library(readr)


portfolio <- 
  read_csv("DATABASE_portfolio.csv")[-60,] 
# obs 60 is an average
portfolio[is.na(portfolio)] <- 0

compute_incremental_roi <- function(prices) {
  # Calculate ROI between each successive data point
  roi <- diff(prices) / head(prices, -1)
  return(roi)
}

roi <- 
  apply(portfolio[,2:10], 2, compute_incremental_roi) %>% 
  as.data.frame() %>% 
  cbind(portfolio$Year[2:59]) %>% 
  .[-1,]

# Define shades of grey and linetypes
grey_shades <- c("black",
                 "grey80", "grey50",
                 "grey20","grey60","grey30")
line_types <- c("solid", "longdash",
                "dotted", "twodash",
                "dotdash", "F1")

melt(roi[,1:6,10]) %>% 
  ggplot(aes(value, color=variable,
             linetype = variable)) +
  geom_density(size=.5) +
  scale_color_manual(values = grey_shades) +
  scale_linetype_manual(values = line_types) +
  theme_minimal() + xlim(-.5,.5) +
   theme(legend.position = "top",
        text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  labs(title = "Densities of Asset Price ROIs",
       x = "Value", y = "Density")


```



```{r  eval=F}
library(tidyverse)
library(kableExtra)
library(bayesAB)


library(fitdistrplus)
fitto <- tibble()
for(i in 1:6){
  fit <- tibble(
    colnames(roi)[i],
    descdist(roi[,i], 
             method='unbiased',
             print=, graph=F)$mean,
     descdist(roi[,i], 
              method='unbiased',
              print=F, graph=F)$sd,
     descdist(roi[,i],
              method='unbiased',
              print=F, graph=F)$skewness,
     descdist(roi[,i],
              method='unbiased',
              print=F, graph=F)$kurtosis
  )
  
  fitto <- rbind(fit, fitto)
}

colnames(fitto) <- c("Asset", 
                     "Mean", "SD", 
                     "Skewness", "Kurtosis")

fitto %>% kable(
  caption = "Fit statistics of
  each asset's time series", 
  digits=3,
  ,booktabs=T
  )




table <- tibble()
   colnames(table) <- c("assetA", 
                        "assetB", 
                        "meanA",
                        "meanB")
for (j in 1:5) {
  for (k in (j + 1):6) {
    roi_test_result <- bayesTest(
      A = roi[, j],
      B = roi[, k],
      priors = c(
        "mu" = 0,
        "lambda" = 1,
        "alpha" = 2,
        "beta" = .05
      ),
      n_samples = 1e+05,
      dist = 'normal'
    )
    
    line <- data.frame(colnames(roi)[j],
                       colnames(roi)[k], mean(t(as.numeric(
      unlist(roi_test_result$posteriors$Mu$A)
    ))), mean(t(as.numeric(
      unlist(roi_test_result$posteriors$Mu$B)
    ))))
    colnames(line) <- c("assetA",
                        "assetB", 
                        "meanA", 
                        "meanB")
    table <- rbind(table, line)
 
  }
}
   
   result <- max(unlist(table[,3:4]))
   winner <- table %>% 
     filter(meanA == result) %>% .[1,c(1,3)]
   winner %>% kable(
     caption="The asset class with the highest ROI", 
     digits=5,
     col.names = 
       c("Asset", "ROI"),
     booktabs=T
       )

```



## R code for Chapter 7


```{r eval=F}

library(tidyverse)
library(bayesAB)

set.seed(123)
A_binom <- rbinom(100, 1, .45)
B_binom <- rbinom(100, 1, .49)

AB1 <- bayesTest(A_binom,
                 B_binom,
                 priors = 
                   c('alpha' = 1, 'beta' = 1),
                 distribution = 'bernoulli')

# summary(AB1)
plot(AB1)

```


```{r eval=F}

set.seed(123)
A_norm <- rnorm(100, 6, 1)
B_norm <- rnorm(100, 5, 5)


AB2 <- bayesTest(
  A_norm,
  B_norm,
  priors = c(
    'mu' = 5,
    'lambda' = 1,
    'alpha' = 3,
    'beta' = 1
  ),
  distribution = 'normal'
)





AB3 <-
  combine(
    AB1,
    AB2,
    f = `*`,
    params = c('Probability', 'Mu'),
    newName = 'Expectation'
  )
##print(AB3)
#summary(AB3)
plot(AB3)




```


```{r eval=F}

set.seed(456)
A_binom <- rbinom(100, 1, .35)
B_binom <- rbinom(100, 1, .30)
A_norm <- rnorm(100, 5.5, 4)
B_norm <- rnorm(100, 5, 3)


## Summarize the posterior means 
#and variances from the first experiment,
#and insert these into the priors of the new A/B test

mu1 <- mean(c(AB1$posteriors$Probability$A,
              AB1$posteriors$Probability$A))
var1 <- var(c(AB1$posteriors$Probability$A,
              AB1$posteriors$Probability$A))

## From Table 2 in Chapter 3 on Priors

alpha1 <- (mu1^2-mu1^3-mu1*var1)/var1
beta1 <-  (mu1-1)*(mu1^2-mu1+var1)/var1


mu2 <- mean(c(AB2$posteriors$Sig_Sq$A,
              AB2$posteriors$Sig_Sq$B))
var2 <- mean(c(AB2$posteriors$Sig_Sq$A,
               AB2$posteriors$Sig_Sq$B))


## From Table 2 in Chapter 3 on Priors
lambda2 <- 1; alpha2 <- 2
mu2 <- mu2
beta2 <- var2


AB1a <- bayesTest(A_binom,
                 B_binom,
                 priors = 
                   c('alpha' = alpha1,
                     'beta' = beta1),
                 distribution = 'bernoulli')

AB2a <- bayesTest(
  A_norm,
  B_norm,
  priors = c(
    'mu' = mu2,
    'lambda' = lambda2,
    'alpha' = alpha2,
    'beta' = beta2
  ),
  distribution = 'normal'
)

AB3 <-
  combine(
    AB1a,
    AB2a,
    f = `*`,
    params = c('Probability', 'Mu'),
    newName = 'Expectation'
  )
#print(AB3)
#summary(AB3)
plot(AB3)



```




```{r eval=F}

library(tidyverse)
library(EnvStats)  # for the Pareto distribution


set.seed(123)
dt <- sample(seq(as.Date('2024-02-01'), 
                 as.Date('2024-05-31'), by = "day"), 10000, replace = TRUE)

NR <- rbernoulli(10000, .3) %>% 
  as.data.frame()
colnames(NR) <- "type"
NR <- NR %>%
  mutate(
    userType = ifelse(type == TRUE,
                      "New Visitor",
                      "Returning Visitor"))
 
NR <- cbind(NR,dt)
NR_T <- NR %>% filter(type == TRUE)
NR_F <- NR %>% filter(type == FALSE)

## Don't differentiate between the userType; 
# let set.seed() preselect whether A or B is the winner
sales_T <- abs(rnorm(nrow(NR_T),
                     25,100)) %>%
  floor()
sales_F <- abs(rnorm(nrow(NR_F),
                     25,100)) %>% 
  floor()

daysSinceLastSession_T <- 0
daysSinceLastSession_F <-
  rpareto(nrow(NR_T), 1, .5) 
daysSinceLastSession_F <- 
  ifelse(
    daysSinceLastSession_F > 100,
    runif(1, 1, 5), 
    daysSinceLastSession_F) %>% floor()

nrt <-   cbind(NR_T, 
               cbind(sales_T, 
                           daysSinceLastSession_T)) %>% 
  rename(
    sales = sales_T,
    daysSinceLastSession = 
      daysSinceLastSession_T)

nrf <-   cbind(NR_F, cbind(sales_F,
                           daysSinceLastSession_F)) %>% 
  rename(
    sales = sales_F,
    daysSinceLastSession =
      daysSinceLastSession_F)

NR <- rbind(nrt,nrf) %>% 
  arrange(dt) %>% 
  rename(date = dt) 


library(bayesAB)


nr_new <- NR %>% 
  filter(userType == "New Visitor")
nr_rtn  <-  NR %>% 
  filter(userType == "Returning Visitor")

AB_new_rtn <- bayesTest(
  nr_new$sales,
  nr_rtn$sales,
  priors = c("shape" = 9, "rate" = 3),
  distribution = 'poisson'
)


plot(AB_new_rtn)


```



```{r eval=F}

library(bayesAB)

nr_short <- NR %>% 
  filter(daysSinceLastSession <= 2)
nr_long <- NR %>% 
  filter(daysSinceLastSession > 2)



AB_delay <- bayesTest(
  nr_short$sales,
  nr_long$sales,
  priors = c("shape" = 9, "rate" = 3),  
  
  distribution = 'poisson'
)

plot(AB_delay)
```


```{r eval=F}


library(bayesAB)

nr_new <- NR %>% filter(userType ==
                          "New Visitor" &
                          daysSinceLastSession <= 2)
nr_rtn <-   NR %>% 
  filter(daysSinceLastSession > 2)
nr_short <- NR %>% 
  filter(daysSinceLastSession <= 2)
nr_long <- NR %>% 
  filter(daysSinceLastSession > 2)

AB_1st <- bayesTest(
  nr_new$sales,
  nr_rtn$sales,
  priors = c("shape" = 9, "rate" = 3),
  ## gives a distribution for number 
  #of days to return centered on 3
  distribution = 'poisson'
)


AB_2nd <- bayesTest(
  nr_short$sales,
  nr_long$sales,
  priors = c("shape" = 9, "rate" = 3),   
  
  distribution = 'poisson'
)



AB_combine <-
  combine(
    AB_1st,
    AB_2nd,
    f = `+`,
    params = 
      c('Lambda', 'Lambda'), 
    newName = 
      'Combined Sales Posterior'
  )


plot(AB_combine)


```


## R code for Chapter 8



```{r eval=F}
# Load necessary libraries
library(ggplot2)
library(gridExtra)
library(knitr)
library(kableExtra)

# Define a function to calculate 
#log-normal parameters, 
#plot PDF and CDF, and compute NPV
plot_lognormal_npv <- 
  function(
    mean_life_expectancy,    
    sd_life_expectancy, 
    color, 
    discount_rate = 0.05, 
    prior_mean, prior_sd) {
  # Convert mean and sd to log-normal parameters
  mu <- 
    log(mean_life_expectancy^2 / 
          sqrt(sd_life_expectancy^2 + 
                 mean_life_expectancy^2))
  sigma <- 
    sqrt(log(1 + (sd_life_expectancy^2 / 
                    mean_life_expectancy^2)))

  # Generate a sequence of ages
  ages <- seq(0, 120, by = 0.1)

  # Calculate the probability density
  #function (PDF) and cumulative density function (CDF)
  pdf_values <- dlnorm(ages, 
                       meanlog = mu, 
                       sdlog = sigma)
  cdf_values <- plnorm(ages,
                       meanlog = mu, 
                       sdlog = sigma)

  # Calculate NPV of $1,000,000 payout 
  #at the time of death
  npv_values <- 10000 * 
    exp(-discount_rate * ages) * pdf_values
  npv <- sum(npv_values) * 
    (ages[2] - ages[1])  
  # Integrate over the age range

  # Bayesian analysis manually
  data <- rlnorm(length(ages), 
                 meanlog = mu, 
                 sdlog = sigma)
  n <- length(data)
  data_mean <- mean(data)
  data_var <- var(data)
  
  # Calculate posterior mean and variance
  post_var <- 1 / 
    (1 / prior_sd^2 + n / data_var)
  post_mean <- post_var * 
    (prior_mean / prior_sd^2 + n * 
       data_mean / data_var)
  posterior_sd <- sqrt(post_var)

  # Calculate NPV from posterior values
  posterior_npv_values <- 
    10000 * exp(-discount_rate * ages) * 
    dlnorm(ages,                             
           meanlog = log(post_mean), 
           sdlog = posterior_sd)
  posterior_npv <- 
    sum(posterior_npv_values) * 
    (ages[2] - ages[1])

  # Create a data frame for plotting
  plot_data <- data.frame(
    Age = ages, 
    PDF = pdf_values, 
    CDF = cdf_values)

  # Plot the PDF
  pdf_plot <- ggplot(plot_data, 
                     aes(x = Age, y = PDF)) +
    geom_line(color = color) +
    ggtitle(paste("Age at Death")) +
    xlim(0,25) + 
    xlab("Age") +
    ylab("Probability Density") +
    theme_minimal()

  # Plot the CDF
  cdf_plot <- ggplot(plot_data, aes(
    x = Age, 
    y = CDF)) +
    geom_line(color = color) +
    ggtitle(paste("Age at Death")) +
    xlab("Age") +
    ylab("Cumulative Probability") +
       xlim(0,25) + 
    theme_minimal()

  list(pdf_plot = pdf_plot, 
       cdf_plot = cdf_plot, 
       npv = npv, 
       posterior_npv = posterior_npv)
}
```



```{r eval=F}


# Scenario 2: Average age of 
#death is 15 with standard deviation 10
scenario2 <- 
  plot_lognormal_npv(mean_life_expectancy = 15, 
                     sd_life_expectancy = 10, 
                     color = "black", 
                     prior_mean = 20,
                     prior_sd = 10)

# Create a data frame for the results
results <- data.frame(
  Scenario = c("Dried Food", "Fresh Food"),
#  Mean_Age = c(5, 15),
#  SD_Age = c(10, 10),
  NPV = 
  round(c(scenario1$npv, 
          scenario2$npv), 2),
  Posterior_NPV = 
  round(c(scenario1$posterior_npv,
          scenario2$posterior_npv), 2)
)



# Display the results in a table
kable(
  results,
  caption = 
    "Comparative Actuarial 
  Costs of Pet Insurance",
  col.names = c(
    "Diet",
    #   "Mean Age",
    #   "SD Age",
    "Frequentist NP Cost of 
    Insurance Claim ($)",
    "Bayesian Posterior NP 
    Cost of Insurance Claim ($)"
  )
) %>%
  kable_styling(full_width = F,
                position = "center"
                ) %>%
  column_spec(1, width = "3cm") %>%
  #  column_spec(2, width = "2cm") %>%
  #  column_spec(3, width = "2cm") %>%
  column_spec(2, width = "3cm") %>%
  column_spec(3, width = "3cm")

# Arrange the plots


library(grid)
library(gridExtra)



grid.arrange(
  scenario1$pdf_plot,
  scenario2$pdf_plot,
  ncol = 2,
  top = textGrob(
    "Mortality Distributions for Dried Food 
    (left) and Fresh Food (right)",
    gp = gpar(fontsize = 14, font = 3)
  )
)




grid.arrange(
  scenario1$cdf_plot,
  scenario2$cdf_plot,
  ncol = 2,
  top = textGrob(
    "Mortality Distributions for Dried Food
    (left) and Fresh Food (right)",
    gp = gpar(fontsize = 14, font = 3)
  )
)



```



```{r  eval=F}
library(tidyverse)
library(readr)
library(bayesAB)

set.seed(123)


girls_names <- read_csv(("girls_names.csv"),
                        col_names = FALSE)[, 2]


sdlog <- 2
meanlog <- 

prof_A <- rpois(1000, 6) # mean 5.5
prof_B <- rpois(1000, 5) # mean 3.7
ron_rate <- rpois(1000, 7) # mean 6.0

prof_A <-ifelse(prof_A > 11, rpois(1000, 6), prof_A)
prof_B <-ifelse(prof_B > 1, rpois(1000, 5), prof_B)
ron_rate <-ifelse(ron_rate > 11, rpois(1000, 7), ron_rate)
prof_A <-ifelse(prof_A > 11, 10, prof_A)
prof_B <-ifelse(prof_B > 11,10, prof_B)
ron_rate <-ifelse(ron_rate > 11, 10, ron_rate)

## this tibble is the "ground truth" 
#dating pool on MatchRate; both sides 
#will in exceptional circumstances assign a rating higher than 10

dating_pool <- tibble(girls_names, prof_A, prof_B, ron_rate)

colnames(dating_pool) <- 
  c("lady", "Profile_A", 
    "Profile_B", 
    "Ladies_rating")
write.csv(dating_pool, 
          file="DATABASE_dating_pool.csv")


# Melt the data for easier plotting with ggplot2
library(reshape2)
dating_pool <- melt(dating_pool)

# Plotting
ggplot(dating_pool, aes(x = value)) +
  geom_histogram(binwidth = 1, fill = "white", 
                 color = "black") +  
  # Adjust binwidth as needed
  facet_wrap(~ variable, scales = "free") +  
  # Separate plot for each variable
  theme_minimal() +
  labs(x = "Ratings (higher is better)", y = "Count") +
  theme(legend.title = element_blank())
```


```{r eval=F}

dating_pool <- 
  read.csv("DATABASE_dating_pool.csv")

# Bayesian A/B test
  bayes_test <- bayesAB::bayesTest(
    prof_A,
    prof_B,
    priors = c("shape" = 1, "rate" = 1),
    distribution = 'poisson'
  )

dating_pool <- data.frame(
  dating_pool$lady,
  bayes_test$posteriors$Lambda$A, 
  bayes_test$posteriors$Lambda$B, 
  dating_pool$Ladies_rating
)
colnames(dating_pool) <- 
  c("lady", 
    "Profile_A", 
    "Profile_B", 
    "Ladies_rating") 

library(reshape2)
dating_pool <- melt(dating_pool)

# Plotting
ggplot(dating_pool, aes(x = value)) +
  geom_histogram(binwidth = .05, 
                 fill = "white", color = "black") +   
  # Adjust binwidth as needed
  facet_wrap(~ variable, scales = "free") + 
  # Separate plot for each variable
  theme_minimal() +
  labs(x = 
         "Ratings (higher is better)", y = "Count") +
  theme(legend.title = element_blank())

  
```




```{r eval=F}


dating_pool <- 
  read.csv("DATABASE_dating_pool.csv")


# Required libraries
library(tidyverse)

set.seed(123)

ron_benchmark <- 8  
# the minimum rating for which 
# either Ron or the ladies will consider dating
lady_benchmark <- 5


max_attempts <- 1000
nights <- tibble()

for (i in 1:1000) {
  i <- j <- 1
  
  # Profile A
  repeat {
    tonight_pool <- 
      dating_pool[sample(nrow(dating_pool),
                         30, 
                         replace = TRUE), ]
    match_A <- subset(tonight_pool, 
                      Profile_A > ron_benchmark & 
                        Ladies_rating > lady_benchmark)
    if (nrow(match_A) > 0) {
      break
    }
    i <- i + 1
    if (i >= max_attempts) {
      print("Max attempts reached
            without finding a 
            match for Profile A.")
      break
    }
  }
  
  # Profile B
  repeat {
    tonight_pool <- 
      dating_pool[sample(nrow(dating_pool), 
                         30, 
                         replace = TRUE), ]
    match_B <- subset(tonight_pool,
                      Profile_B > ron_benchmark & 
                        Ladies_rating > lady_benchmark)
    if (nrow(match_B) > 0) {
      break
    }
    j <- j + 1
    if (j >= max_attempts) {
      print("Max attempts 
            reached without finding a match for Profile B.")
      break
    }
  }
  
  nights <- rbind(nights, tibble(Prof_A = i, Prof_B = j))
}

colnames(nights) <- c("Profile A", "Profile B")

# Reshaping and plotting
nights_long <- pivot_longer(nights, everything(),
                            names_to = "profile", 
                            values_to = "nights")
ggplot(nights_long, aes(x = nights)) +
  geom_histogram(binwidth = 1, 
                 fill = "white", color = "black") +
  facet_wrap(~profile, scales = "free") +
  theme_minimal() +
  labs(x = "Number of Nights for
       Ron to Meet his Match", 
       y = "How Often (Number of
       Nights out of 1000)") +
  theme(legend.title = element_blank()) +
  scale_y_continuous(trans = 'log10')




```




## R code for Chapter 9


```{r eval=F}
library(wordcloud)
library(tidyverse)
library(tidytext)
library(tokenizers)
library(tm)
library(readr)
library(knitr)

reviews <- read_csv("DATABASE_reviews.csv") %>%
  select(UserId, Score, Text)  

set.seed(123)
mkt <- sample(seq_len(nrow(reviews)), 
              size = nrow(reviews) / 2)
A <- reviews %>% filter(Score < 4) 
B <- reviews %>% filter(Score >= 4) 

# write.csv(A, B, reviews, 'NLP_reviews.csv')
  
```

NLP Analysis

```{r eval=F}
# Load necessary libraries
library(dplyr)
library(tidyr)
library(tokenizers)
library(tidytext)

# Check and preprocess text data
preprocess_text <- function(data) {
  if (is.data.frame(data) && "Text" %in% names(data)) {
    data$Text <- as.character(data$Text)
    Encoding(data$Text) <- "UTF-8"
    text_stuff <- tokenize_words(data$Text) %>%
      unlist() %>%
      as.data.frame(stringsAsFactors = FALSE)
    colnames(text_stuff) <- "word"
    return(text_stuff)
  } else {
    stop("Input data must be a
         data frame with a 'Text' column")
  }
}

# Compute net sentiment
compute_net_sentiment <- function(text_stuff) {
  sentiments <- get_sentiments("nrc")
  net_sentiment <- text_stuff %>%
    inner_join(sentiments, by = "word") %>%
    count(sentiment) %>%
    pivot_wider(names_from = sentiment, 
                values_from = n, 
                values_fill = list(n = 0)) %>%
    mutate(net_positive = positive - negative,
           proportion_positive = 
             positive / negative - 1)
  return(net_sentiment)
}



# Process data A and B if they are correctly formatted
text_stuff_A <- preprocess_text(A)
net_sentiment_A <- compute_net_sentiment(text_stuff_A)

text_stuff_B <- preprocess_text(B)
net_sentiment_B <- compute_net_sentiment(text_stuff_B)

# Combine results and create a summary table
net_sentiment <- rbind(net_sentiment_A, net_sentiment_B)
product <- data.frame(Product = c("A", "B"),
                      stringsAsFactors = FALSE)
net_sentiment <- cbind(product, net_sentiment)

# Print the results
print(net_sentiment)



```



```{r eval=F}
# Load necessary libraries
library(dplyr)
library(tidyr)
library(tokenizers)
library(tidytext)

# Check and preprocess text data
preprocess_text <- function(data) {
  if (is.data.frame(data) && "Text" %in% names(data)) {
    Encoding(data$Text) <- "UTF-8"
    text_stuff <- tokenize_words(data$Text) %>%
      unlist() %>%
      as.data.frame()
    colnames(text_stuff) <- "word"
    return(text_stuff)
  } else {
    stop("Input data must be a data 
         frame with a 'Text' column")
  }
}

# Compute net sentiment
compute_net_sentiment <- function(text_stuff) {
  net_sentiment <- text_stuff %>%
    inner_join(get_sentiments("nrc"), by = "word") %>%
    count(sentiment) %>%
    spread(sentiment, n, fill = 0) %>%
    mutate(net_positive = positive - negative,
           proportion_positive = 
             positive / negative - 1)
  return(net_sentiment)
}

# Process data A and B if they are correctly formatted
text_stuff_A <- preprocess_text(A)
net_sentiment_A <- 
  compute_net_sentiment(text_stuff_A)

text_stuff_B <- preprocess_text(B)
net_sentiment_B <- 
  compute_net_sentiment(text_stuff_B)

# Combine results and create a summary table
net_sentiment <- 
  rbind(net_sentiment_A, net_sentiment_B)
product <- data.frame(Product = c("A", "B"))
net_sentiment <- 
  cbind(product, net_sentiment)

# Print the results
print(net_sentiment)




```



```{r eval=F}

# Load necessary libraries
library(dplyr)
library(tidyr)
library(tokenizers)
library(tidytext)

# Function to check and preprocess text data
preprocess_text <- function(data) {
  # Validate that the input is a data frame
  if (!is.data.frame(data)) {
    stop("Input data must be a data frame.")
  }
  
  # Validate that the data frame contains a 'Text' column
  if (!"Text" %in% names(data)) {
    stop("Data frame must contain a 'Text' column.")
  }
  
  # Set the text encoding to UTF-8 and preprocess
  Encoding(data$Text) <- "UTF-8"
  text_stuff <- tokenize_words(data$Text) %>%
    unlist() %>%
    as.data.frame()
  colnames(text_stuff) <- "word"
  return(text_stuff)
}

# Function to compute net sentiment
compute_net_sentiment <- function(text_stuff) {
  net_sentiment <- text_stuff %>%
    inner_join(get_sentiments("nrc"), by = "word") %>%
    count(sentiment) %>%
    spread(sentiment, n, fill = 0) %>%
    mutate(net_positive = positive - negative,
           proportion_positive = positive / negative - 1)
  return(net_sentiment)
}


# Process text data from A, ensuring it's properly formatted
if (is.data.frame(A) && "Text" %in% names(A)) {
  text_stuff_A <- preprocess_text(A)
  net_sentiment_A <- 
    compute_net_sentiment(text_stuff_A)
  
  # Combine results and create a summary table
  net_sentiment <- rbind(net_sentiment_A) 
  # Add more datasets as needed
  product <- data.frame(Product = c("A")) 
  # Extend with more products as needed
  net_sentiment <- cbind(product, net_sentiment)

  # Print the results
  print(net_sentiment)
} else {
  print("Error: The data 'A' is not a 
        data frame or lacks a 'Text' column.")
}


```


WordClouds

```{r eval=F}


A <- A[ceiling(runif(10000,0,nrow(A))),]  
B <- B[ceiling(runif(10000,0,nrow(B))),]

cat("Wordcloud for Product A")

Encoding(A$Text) = "UTF-8" 
text_stuff <- tokenize_words(A$Text) %>% 
  unlist() %>%
  as.data.frame()

colnames(text_stuff) <- "word"

stuff_sentiment <-
  text_stuff %>% 
  inner_join(get_sentiments("nrc"), by = "word")

stp <- get_stopwords()
stp <- rbind(stp, "good", "food", "love", 
             "sweet", "chocolate", "bad", "sugar",
             "found", "weight", "smell", "money")

frequency_A <- text_stuff %>% anti_join(stp) %>% 
  inner_join(stuff_sentiment) %>% count(word)
frequency_A  %>%
  with(wordcloud(
    word,
 colors = "black",
    rot.per = 0.15,
    n,
    max.words = 70
  ))

cat("Wordcloud for Product B")

Encoding(A$Text) = "UTF-8" 
text_stuff <- tokenize_words(B$Text) %>% 
  unlist() %>%
  as.data.frame()

colnames(text_stuff) <- "word"

stuff_sentiment <-
  text_stuff %>%
  inner_join(get_sentiments("nrc"), 
             by = "word")

frequency_B <- text_stuff %>% anti_join(stp) %>% 
  inner_join(stuff_sentiment) %>% count(word)
frequency_B  %>%
  with(wordcloud(
    word,
    colors = "black",
    rot.per = 0.15,
    n,
    max.words = 70
  ))



```


Top 10 Words for Each Sentiment    

```{r eval=F}


cat("Top 10 Words by Sentiment for Product A")

Encoding(A$Text) = "UTF-8" 
text_stuff <- tokenize_words(A$Text) %>% 
  unlist() %>%
  as.data.frame()

colnames(text_stuff) <- "word"


stuff_sentiment <-
text_stuff %>%  anti_join(stp) %>% 
  inner_join(get_sentiments("nrc"), by = "word")

text_counts <- stuff_sentiment %>%
count(word, sentiment, sort = TRUE) %>%
ungroup()

text_counts %>%
  group_by(sentiment) %>%
  top_n(10) %>%
  ungroup() %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(word, n, fill = sentiment)) +
  scale_fill_grey(start = 0, end = .2) +
  theme_minimal() +  
  # Optional: starts with a minimal theme
  theme(axis.text.
        x = element_text(angle = 45, hjust = 1)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ sentiment, scales = "free_y") +
  labs(y = "Contribution to sentiment", x = NULL) +
  coord_flip()


cat("Top 10 Words by Sentiment for Product B")

Encoding(B$Text) = "UTF-8" 
text_stuff <- tokenize_words(B$Text) %>% 
  unlist() %>%
  as.data.frame()

colnames(text_stuff) <- "word"


stuff_sentiment <-
text_stuff %>% anti_join(stp) %>%  
  inner_join(get_sentiments("nrc"), by = "word")

text_counts <- stuff_sentiment %>%
count(word, sentiment, sort = TRUE) %>%
ungroup()

text_counts %>%
  group_by(sentiment) %>%
  top_n(10) %>%
  ungroup() %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(word, n, fill = sentiment)) +
  scale_fill_grey(start = 0, end = .2) +
  theme_minimal() +  
  # Optional: starts with a minimal theme
  theme(axis.text.x = element_text(
    angle = 45, hjust = 1)
    ) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ sentiment, scales = "free_y") +
  labs(y = "Contribution to sentiment", x = NULL) +
  coord_flip()


```



```{r eval=F}


library(tidyverse)
library(tidytext)
library(tokenizers)
library(tm)
library(readr)
library(knitr)


Encoding(A$Text) = "UTF-8"
text_stuff <- tokenize_words(A$Text) %>%
  unlist() %>%
  as.data.frame()
colnames(text_stuff) <- "word"
stuff_sentiment_A <-
  text_stuff %>%  
  inner_join(get_sentiments("nrc"), by = "word")  %>% 
  count(sentiment, word)

sum_sent <- stuff_sentiment_A %>%  
  group_by(sentiment) %>% 
  summarize(mu = mean(n), sig = sd(n))
sent_stats_A <-
  data.frame(sum_sent$sentiment,
             sum_sent$mu , 
             sum_sent$sig)
sent_stats_A


Encoding(B$Text) = "UTF-8"
text_stuff <- tokenize_words(B$Text) %>%
  unlist() %>%
  as.data.frame()
colnames(text_stuff) <- "word"
stuff_sentiment_B <-
  text_stuff %>% 
  inner_join(get_sentiments("nrc"), 
             by = "word")  %>% 
  count(sentiment, word)

sum_sent <- stuff_sentiment_B %>%  
  group_by(sentiment) %>% 
  summarize(mu = mean(n),
            sig = sd(n))
sent_stats_B <- 
  data.frame(sum_sent$sentiment,
             sum_sent$mu , 
             sum_sent$sig)
sent_stats_B

```


```{r eval=F}

library(bayesAB)

snt <- unique(stuff_sentiment_A$sentiment)
sent_stats <- data.frame()
for (i in snt) {
  sent_A <- stuff_sentiment_A %>% 
    filter(sentiment == i)
  sent_B <- stuff_sentiment_B %>%
    filter(sentiment == i)
  AB1 <- bayesTest(
    sent_A$n,
    sent_B$n,
    priors = c(
      "mu" = 0,
      "lambda" = 1,
      "alpha" = 1,
      "beta" = 1
    ),
    distribution = 'normal'
  )
cat("Sentiment = ", i)  
#print(AB1)
#summary(AB1)
#print(plot(AB1))

# Load the necessary library
library(ggplot2)

# Example data vectors
Choice_A <- AB1$posteriors$Mu$A  
# Generate some normal data
Choice_B <- AB1$posteriors$Mu$B  
# Generate some normal data 

# Create a data frame for ggplot
data <- data.frame(value = c(data1, data2),
                   group = 
                     factor(rep(c("Choice A",
                                  "Choice B"), 
                                each = 100)))

# Create the density plot
p <- ggplot(data, aes(x = value, group = group)) +
  geom_density(aes(linetype = group), size = 1) +  
  # Use different linetypes for each group
  scale_linetype_manual(values = c("solid", "dashed")) +  
  # Assign solid and dashed lines
  facet_wrap(~ group, scales = "free_y") + 
  # Facet by group, allowing free scales on the y-axis
  labs(title = "Density Plot by Group",
       x = i, y = "Density") +
  theme_minimal() +  
  # Using a minimal theme for better visibility
  theme(text = element_text(color = "black"),  
        # Ensure text is black for contrast
        panel.grid.major = element_blank(),   
        # Remove major grid lines
        panel.grid.minor = element_blank(),  
        # Remove minor grid lines
        panel.background = element_blank(),   
        # Remove panel background
        plot.background = element_blank(),    
        # Transparent plot background
        strip.background = element_blank(),  
        # Transparent facet strip background
        strip.text = element_text(color = "black")) 

print(p)


}

```




```{r eval=F}

# Load the necessary library
library(ggplot2)

# Example data vectors
data1 <- rnorm(100, mean = 50, sd = 10) 
# Generate some normal data
data2 <- rnorm(100, mean = 60, sd = 15)  
# Generate some normal data with different parameters

# Create a data frame for ggplot
data <- data.frame(value = c(data1, data2), group =
                     factor(rep(c("Data1", 
                                  "Data2"), 
                                each = 100)))

# Create the density plot
ggplot(data) +
  geom_density(aes(x = value, fill = group, 
                   color = group), alpha = 0.5) +
  scale_fill_manual(values = c("blue", "red")) +  
  # Set custom fill colors for the groups
  scale_color_manual(values = c("blue", "red")) + 
  # Set custom line colors for the groups
  labs(title = "Density Plot",
       x = "Value",
       y = "Density") +
  theme_minimal()  
# Using a minimal theme for better visibility


```



```{r eval=F}

# Load the necessary library
library(ggplot2)

# Example data vectors
Choice_A <- AB1$posteriors$Mu$A 
# Generate some normal data
Choice_B <- AB1$posteriors$Mu$B 
# Generate some normal data with different parameters

# Create a data frame for ggplot
data <- data.frame(value = c(data1, data2),
                   group =
                     factor(rep(c("Choice A", "Choice B"), 
                                each = 100)))

# Create the density plot
ggplot(data, aes(x = value, group = group)) +
  geom_density(aes(linetype = group), size = .3) +  
  # Different linetypes for each group
  scale_linetype_manual(values =
                          c("solid", "dashed")) +  
  # Assign solid and dashed lines
  labs(title = "Density Plot",
       x = "Value", 
       y = "Density") +
  theme_minimal() +  
  # Using a minimal theme for better visibility
  theme(text = element_text(color = "black"),  
        # Ensure text is black for contrast
        panel.grid.major = element_blank(),  
        # Remove major grid lines
        panel.grid.minor = element_blank(),  
        # Remove minor grid lines
        panel.background = element_blank(),  
        # Remove panel background
        plot.background = element_blank())    
# Transparent plot background


```








```{r eval=F}

library(bayesAB)
library(readr)

priors <- 
  read_csv(
    "DATABASE_priors_for_NLP_case_study.csv"
    )[1:10,c(2,5:8)]



snt <- unique(stuff_sentiment_A$sentiment)
sent_stats <- data.frame()
for (i in c(
  "anger",
  joy,
  "anticipation",
  "disgust",
  "fear", 
  "sadness",
  "surprise",
  "trust",
  "negative",
  "positive"
  )) {
  sent_A <-
    stuff_sentiment_A %>% filter(sentiment == i)
  sent_B <- 
    stuff_sentiment_B %>% filter(sentiment == i)
  prr <- priors %>% filter(Sentiment == i)
  AB1 <- bayesTest(
    sent_A$n,
    sent_B$n,
    priors = c(
      "mu" = prr$mu,
      "lambda" = prr$lambda,
      "alpha" = prr$alpha,
      "beta" = prr$beta
    ),
    distribution = 'normal'
  )
  

cat("Sentiment = ", i)  
#print(AB1)
#summary(AB1)
#print(plot(AB1))

# Load the necessary library
library(ggplot2)

# Example data vectors
data1 <- AB1$posteriors$Mu$A 
# Generate some normal data
data2 <- AB1$posteriors$Mu$A  
# Generate some normal data with different parameters

# Create a data frame for ggplot
data <- data.frame(value = c(data1, data2),
                   group = 
                     factor(rep(c("Data1",
                                  "Data2"), 
                                each = 100)))

# Create the density plot
ggplot(data, aes(x = value, 
                 fill = group, 
                 alpha = 0.5)) +
  geom_density(adjust = 1.5) +  
  # Adjust parameter on how smooth you want the curve to be
  scale_fill_manual(values = c("blue", "red")) +  
  # Set custom colors for the groups
  scale_alpha_manual(values = c(0.5, 0.5)) + 
  # Set transparency level
  labs(title = "Density Plot", 
       x = "Value", 
       y = "Density") +
  theme_minimal()  
# Using a minimal theme for better visibility


}

```


```{r eval=F}

library(kableExtra)
library(tidyverse)
library(tidytext)
library(tokenizers)
library(tm)
library(readr)
library(knitr)


reviews <- read_csv("DATABASE_reviews.csv") %>% 
  select(UserId, Score, Text)

set.seed(123)
mkt <- sample(seq_len(nrow(reviews)), 
              size = nrow(reviews) / 2)
A <- reviews %>% filter(Score < 4)
B <- reviews %>% filter(Score >= 4)

Encoding(A$Text) = "UTF-8"
text_stuff <- tokenize_words(A$Text) %>%
  unlist() %>%
  as.data.frame()
colnames(text_stuff) <- "word"
stuff_sentiment_A <-
  text_stuff %>%  
  inner_join(get_sentiments("nrc"), by = "word")  %>% 
  count(sentiment, word)

sum_sent <- 
  stuff_sentiment_A %>% 
  group_by(sentiment) %>% 
  summarize(mu = mean(n), sig = sd(n))
sent_stats_A <-
  data.frame(sum_sent$sentiment,
             sum_sent$mu ,
             sum_sent$sig)
sent_stats_A %>% kbl(
  caption="Sentiment A Statistics")


Encoding(B$Text) = "UTF-8"
text_stuff <- tokenize_words(B$Text) %>%
  unlist() %>%
  as.data.frame()
colnames(text_stuff) <- "word"
stuff_sentiment_B <-
  text_stuff %>%
  inner_join(get_sentiments("nrc"), by = "word")  %>% 
  count(sentiment, word)

sum_sent <- stuff_sentiment_B %>% 
  group_by(sentiment) %>% 
  summarize(mu = mean(n), sig = sd(n))
sent_stats_B <- 
  data.frame(sum_sent$sentiment, 
             sum_sent$mu , 
             sum_sent$sig)
sent_stats_B %>% kbl(caption="Sentiment B Statistics")

```







## R code for Chapter 10


```{r eval=F}

library(bayesAB)
library(tidyverse)

bs <- read_csv("kildare_bs.csv")


pre_theft <-  bs %>% 
  filter(Year < 2015) %>% 
  select(Equity) %>% 
  unlist() %>% 
  as.integer()
post_theft <- bs %>% 
  filter(Year >= 2015) %>% 
  select(Equity) %>% 
  unlist() %>% 
  as.integer()

AB_NPV <- bayesTest(
    pre_theft,
    post_theft,
    priors = c("shape" = 10, 
               "rate" = 10/mean(bs$Equity)),  
    distribution = 'poisson'
  )

plot(AB_NPV)


```

## R code for Chapter 11


```{r eval=F}

library(readr)
library(bayesAB)
library(tidyverse)



fin <- read_csv("kildare_inc_stmt.csv")
deflator <- .04  # assume 4% annual inflation

efin <- dfin <- fin %>%
  mutate(
    deflate = (1 + deflator)^(Year-2023),
    Revenue = trunc(Revenue / deflate),
    CGS = trunc(CGS / deflate),
    GP = trunc(GP / deflate),
    RD = trunc(RD / deflate),
    OE = trunc(OE / deflate),
    NI = trunc(NI / deflate)
  )


pre_theft <-  efin %>% 
  filter(Year < 2015) %>% 
  select(NI) %>% 
  unlist() %>% 
  as.numeric()
post_theft <- efin %>% 
  filter(Year >= 2015) %>% 
  select(NI) %>% 
  unlist() %>%
  as.numeric()

AB_profit <- bayesTest(
    pre_theft,
    post_theft,
    priors = c("shape" = 10, "rate" = 10/mean(efin$NI)),  
    distribution = 'poisson'
  )

plot(AB_profit)


```

Bandit with Bernoulli distributed data and Thompson sampling


```{r eval=F}


library(ggplot2)
library(gridExtra)

# Function to simulate bandit arms 
# with Bernoulli rewards

make_bandits <- function(params) {
  pull <- function(arm) {
    function() {
      rbinom(1, 1, params[arm + 1])
    }
  }
  
  list(pull = pull, 
       num_bandits = length(params))
}

# Bayesian strategy for bandit selection
bayesian_strategy <- 
  function(pull, num_bandits) {
  num_rewards <- numeric(num_bandits)
  num_trials <- numeric(num_bandits)
  
  list(
    play = function() {
      choice <- 
        which.max(rbeta(
          num_bandits, 
          2 + num_rewards, 
          2 + num_trials - num_rewards)) - 1
      reward <- pull(choice)()
      num_rewards[choice + 1] <- 
        num_rewards[choice + 1] + reward
      num_trials[choice + 1] <- 
        num_trials[choice + 1] + 1
      list(choice = choice, 
           reward = reward, 
           num_rewards = num_rewards, 
           num_trials = num_trials)
    }
  )
}

# Plot posterior distributions
plot_posteriors <- function(num_rewards, 
                            num_trials, 
                            n, 
                            title_suffix) {
  x <- seq(0, 1, length.out = 100)
  data <- expand.grid(x = x, 
                      Bandit = seq_along(num_rewards))
  densities <- mapply(
    function(a, b) dbeta(x, a + 2, b - a + 2),
    a = num_rewards,
    b = num_trials,
    SIMPLIFY = FALSE
  )
  
  # Flatten the list 
  # and attach it to the dataframe
  data$y <- unlist(densities)
  
  ggplot(data, aes(x = x, 
                   y = y, 
                   group = Bandit)) +
    geom_line(aes(linetype = as.factor(Bandit)), 
              color = "grey20") +
    geom_area(aes(fill = as.factor(Bandit)), 
              alpha = 0.3) +
    scale_linetype_manual(
      values = c("solid", "dashed", "dotted")) +
    scale_fill_manual(
      values = c("grey10", "grey30", "grey50")) +
    labs(title = paste('Posterior after', 
                       n, 
                       'pulls',
                       title_suffix),
         x = 'Probability', 
         y = 'Density') +
    theme_minimal() +
    theme(legend.title = element_blank())
}

# Setup and play the game
bandit_setup <- 
  make_bandits(c(0.2, 0.5, 0.7))
num_bandits <- 
  bandit_setup$num_bandits
play_function <- 
  bayesian_strategy(bandit_setup$pull, 
                    num_bandits)$play

pulls_to_plot <- c(1, 2, 10, 20, 50, 100)

# Plot at specific points
plots <- lapply(pulls_to_plot, function(n) {
  for (i in 1:n) {
    result <- play_function()
  }
  plot_posteriors(result$num_rewards, 
                  result$num_trials, 
                  n, 
                  ifelse(n > 1, 'pulls', 'pull'))
})

# Display plots
gridExtra::grid.arrange(grobs = plots, ncol = 2)




```

Bandit code with Count Distributions


```{r  eval=F}

library(ggplot2)
library(reshape2)  # For melting data frames
library(gridExtra) # For arranging plots

set.seed(123)

# Create bandits with Poisson-distributed rewards
make_bandits <- function(params) {
  pull <- function(arm) {
    function() {
      rpois(1, lambda = params[arm + 1])
    }
  }
  
  list(pull = pull, num_bandits = length(params))
}

# Bayesian strategy, Gamma priors
bayesian_strategy <- function(pull, num_bandits) {
  num_rewards <- rep(1, num_bandits)
  num_trials <- rep(1, num_bandits)
  
  list(
    play = function() {
      choice <- 
        which.max(rgamma(num_bandits, 
                                 shape = num_rewards, 
                                 rate = num_trials)) - 1
      reward <- pull(choice)()
      num_rewards[choice + 1] <-
        num_rewards[choice + 1] + reward
      num_trials[choice + 1] <-
        num_trials[choice + 1] + 1
      list(
        choice = choice,
        reward = reward,
        num_rewards = num_rewards,
        num_trials = num_trials
      )
    }
  )
}

# Function to plot the posterior 
plot_posteriors <- 
  function(num_rewards, 
           num_trials, 
           n) {
  x <- seq(0, 10, length.out = 100)
  plot_data <- data.frame(x = x)
  
  # Compute densities
  for (i in seq_along(num_rewards)) {
    plot_data[[paste0("Bandit", i)]] <-
      dgamma(x, shape = num_rewards[i],
             rate = num_trials[i])
  }
  
  # Reshape for plotting
  plot_data_melted <- melt(
    plot_data,
    id.vars = "x",
    variable.name = "Bandit",
    value.name = "Density"
  )
  
  # Plot configurations
  linetypes <- c("solid", "dashed", "dotted")
  colors <- c("black", "black", "black")
  
  p <- ggplot(plot_data_melted, aes(x = x, 
                                    y = Density, 
                                    group = Bandit)) +
    geom_line(aes(linetype = Bandit), 
              size = .4, 
              color = colors[as.integer(as.factor(
                plot_data_melted$Bandit))]) +
    geom_area(aes(fill = Bandit), alpha = 0.1) +
    scale_linetype_manual(values = linetypes) +
    scale_fill_manual(values = colors) +
    labs(
      title = sprintf("Posterior 
                      distributions 
                      after %d pulls", 
                      n),
      x = "Reward size",
      y = "Density"
    ) +
    theme_minimal() +
    theme(legend.title = element_blank())
  
  return(p)
}

# Initialize bandits and play function
bandit_setup <- make_bandits(c(4, 4.5, 5))
num_bandits <- bandit_setup$num_bandits
play_function <- bayesian_strategy(
  bandit_setup$pull,
  num_bandits)$play

# Collect plots for specific pulls
pulls_to_plot <- c(1, 2, 10, 20, 50, 100)
plots <- list()

for (i in 1:max(pulls_to_plot)) {
  result <- play_function()
  if (i %in% pulls_to_plot) {
    plots[[length(plots) + 1]] <-
      plot_posteriors(
        result$num_rewards, 
        result$num_trials, i)
  }
}

# Arrange and display all plots in a grid
do.call(grid.arrange, c(plots, ncol = 2))

```

Shiny example for Real-time interactive implementation

```{r  eval=F}

library(shiny)
library(ggplot2)
library(gridExtra)

# Define UI for application 
ui <- fluidPage(
  titlePanel("Bayesian Bandit Simulation"),
  sidebarLayout(
    sidebarPanel(
      helpText("Simulate bandit arms with 
               Bayesian updates on 
               reward probabilities."),
      numericInput("bandit1", 
                   "Probability of 
                   reward for Bandit 1:", 
                   0.2, 
                   min = 0, 
                   max = 1),
      numericInput("bandit2", 
                   "Probability of
                   reward for Bandit 2:", 
                   0.5,
                   min = 0, 
                   max = 1),
      numericInput("bandit3", 
                   "Probability of 
                   reward for Bandit 3:",
                   0.7, 
                   min = 0, 
                   max = 1),
      sliderInput("pulls", "Number of pulls:", 
                  min = 1, 
                  max = 500, 
                  value = 100),
      actionButton("simulate", "Simulate")
    ),
    mainPanel(
      plotOutput("plots")
    )
  )
)

# Define server logic 
server <- function(input, output) {
  observeEvent(input$simulate, {
    isolate({
      params <- c(input$bandit1, 
                  input$bandit2,
                  input$bandit3)
      num_bandits <- length(params)
      num_rewards <- numeric(num_bandits)
      num_trials <- numeric(num_bandits)
      
      pull <- function(arm) {
        function() {
          rbinom(1, 1, params[arm + 1])
        }
      }
      
      for (i in 1:input$pulls) {
        choice <- which.max(rbeta(
          num_bandits, 
          2 + num_rewards, 
          2 + num_trials - num_rewards)) - 1
        reward <- pull(choice)()
        num_rewards[choice + 1] <- 
          num_rewards[choice + 1] + reward
        num_trials[choice + 1] <- 
          num_trials[choice + 1] + 1
      }
      
      plots <- 
        lapply(seq_len(num_bandits), 
               function(bandit) {
        data <- data.frame(
          x = seq(0, 1, length.out = 100),
          y = dbeta(seq(0, 1,
                        length.out = 100), 
                    num_rewards[bandit] + 2,
                    num_trials[bandit] - 
                      num_rewards[bandit] + 2)
        )
        ggplot(data, aes(x = x, y = y)) +
          geom_line() +
          geom_area(fill = "gray", alpha = 0.3) +
          labs(title = sprintf(
            "Posterior for Bandit %d after %d pulls",
            bandit, input$pulls),
               x = 'Probability', y = 'Density') +
          theme_minimal()
      })
      
      output$plots <- renderPlot({
        do.call(gridExtra::grid.arrange, 
                c(plots, ncol = 1))
      })
    })
  })
}

# Run the application 
shinyApp(ui = ui, server = server)
```