Authors@R: person(“Annice”, “Najafi”, email = “annicenajafi27@gmail.com”)
Summer 2024
This is an R implementation of the gambling model with changing win probabilities from Chapter 1, section 4 (“A Gambling Model With Changing Win Probabilities”) of the “Introduction to Stochastic Dynamic Programming” book by Sheldon Ross.
A gambler that has x fortune, makes n bets wins with probability p At each stage the gambler has p probability of winning, with alpha probability of winning for previous stages taken from a distribution F the gambler then decides to bet y from 0 to his current fortune (x) goal is to maximize the expectation of the log of his fortune Note: utility function taken to be u(x) = log(x) like before, adding 1 to prevent errors. NOTE: alpha is not the percentage we are betting but rather the probability of winning at the previous stage, don’t confuse this with the alpha in the gambling model in section 2.
Let’s start by loading relevant libraries
#Load relevant libraries
library(ggplot2)
library(truncnorm)
# Define the probability density function, going with the truncated normal dist
#' probability density function for the truncated normal dist
#'
#' @param alpha the winning probability
#' @return the pdf
#' @export
f <- function(alpha) {
return(dtruncnorm(alpha, a = 0, b = 1, mean = 0.5, sd = 0.1))
}
#Great now lets define the utility function
#' Utility function
#'
#' @param x
#' @return log(x+1), adding 1 because if the gambler bets all of his fortune
#' then we will get an error because log(0) is undefined
#' @export
u <- function(x){
return(log(x+1))
}
#' Recursive function that returns the maximum expected reward
#'
#' @param n how many times we are betting
#' @param x the fortune
#' @param p the probability of winning
#'
#' @return the maximal expected reward
#' @export
V <- function(n, x, p) {
# Obviously first we define the base case
if (n == 0) {
return(u(x))
}
alpha <- rtruncnorm(1, a = 0, b = 1, mean = 0.5, sd = 0.1) # This is your alpha taken randomly
hold.integrand_win <- function(alpha) {
V(n - 1, (x + y), alpha) * f(alpha)
}
hold.integrand_lose <- function(alpha) {
V(n - 1, (x - y), alpha) * f(alpha)
}
max_y <- 0 # Initialize max_y
#Find the best y:
for (y in 0:x) {
val <- p * integrate(hold.integrand_win, 0, 1)$value +
(1 - p) * integrate(hold.integrand_lose, 0, 1)$value
# Use max to find the maximum value
max_y <- max(max_y, val)
}
return(max_y) #Return the maximal reward found through searching all y values
}
# Initial fortune
x_initial <- 10
# Initial win probability
p_initial <- 0.5
# Number of plays
n <- 3
c()->return_vect
for(p_initial in seq(0.5, 1, 0.05)){
# Calculate maximal expected final utility
max_expected_utility <- V(n, x_initial, p_initial)
return_vect <- c(return_vect, max_expected_utility)
}
df <- data.frame(p= seq(0.5, 1, 0.05), return=return_vect)
ggplot(df, aes(x=p, y=return_vect))+geom_line(size=2, color = "#538392")+
geom_point(color="#538392", shape=8, stroke=3)+
labs(x="Probability of winning", y="Return")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="A")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
hold.vect <- c()
x<-10
p<-0.8
n<-3
for (y in seq(0, x, 1)){
hold.val <- V(n - 1, x + y, p) - V(n - 1, x - y, p)
hold.vect <- c(hold.vect, hold.val)
}
y_increasing <- data.frame(y = seq(0, x, 1), equation = hold.vect)
ggplot(y_increasing, aes(x=y, y=equation))+geom_line(size=2, color = "#538392")+
geom_point(color="#538392", shape=8, stroke=3)+
labs(x="y", y="Return")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="B")