library(ggplot2)
options( warn =-1)
############################## configuration ##############################
# set working directory
setwd("C:/Users/robin/Desktop/Monopolist")
# do spline? if false, we'll do line
is_spline <- TRUE
# if do_spline, which method?
#spline_method <- "monoH.FC" # straight lines and bends through points
#spline_method <- "periodic" # bends, but not necessarly through points, flattens toward sine
spline_method <- "fmm" # default, alwyas bends, through points
# if not do_spline, which method?
#nspline_method <- "linear" #lines
nspline_method <- "constant" #blocks
# Number of customers %% 200
customers <- 20000
# LiF settings
amount_features_start <- 20
#TT <- customers
#T <- 1.35
#w <- (2*pi)/T
#A <- 21
#gamma <- 0.04
#window.length <- 1
#TT <- customers
#T <- 8
#w <- (2*pi)/T
#A <- 20
#gamma <- 0.8
#window.length <- 1
TT <- customers
T <- 8
w <- (2*pi)/T
A <- 20
gamma <- 0.9
window.length <- 1
############################## helper functions #########################
push <- function(A, value){
if(any(is.na(A))){
i <- sum(!is.na(A))
A[i+1] <- value
} else {
A <- c(A,value)
A <- A[-1]
}
return(A)
}
slideFunct <- function(data, window, step){
total <- length(data)
spots <- seq(from=1, to=(total-window), by=step)
result <- vector(length = length(spots))
for(i in 1:length(spots)){
result[i] <- mean(data[spots[i]:(spots[i]+window)])
}
return(result)
}
############################## sketch some graphs #####################
x_decoy_features <- c(30,60) # A has 30 features B has 60 features
y_decoy_price <- c(30,60) # A costs 30 B costs 60
plot(x_decoy_features, y_decoy_price, xlim = c(0, 90), ylim = c(90,0), col="red")
text(x_decoy_features+0.2, y_decoy_price-4, labels = c("A","B"), col="red")
abline(h=60,lty = 2)
text(25, 65, labels = "Decoy values of C")
text(25, 71, labels = "along dashed line")
x1 <- c( 0, 30, 60, 90 )
f1 <- c( .4, .4, .2, .2 )
x2 <- c( 0, 30, 60, 90 )
f2 <- c( .4, .4, .7, .0 )
x3 <- c( 0, 30, 60, 90 )
f3 <- c( .2, .2, .1, .8 )
if (is_spline) {
s1 <- splinefun(x1, f1,method=spline_method)
plot(x1, f1, ylim = c(0, 1))
curve(s1(x), add = TRUE, col = 2, n = 1001)
s2 <- splinefun(x2, f2,method=spline_method)
plot(x2, f2, ylim = c(0, 1))
curve(s2(x), add = TRUE, col = 2, n = 1001)
s3 <- function(x){
return( 1 - s1(x) - s2(x) )
}
plot(x3, f3, ylim = c(0, 1))
curve(s3(x), add = TRUE, col = 2, n = 1001)
} else {
s1 <- approxfun(x1, f1,method = nspline_method)
plot(x1, f1, ylim = c(0, 1))
curve(s1(x), add = TRUE, col = 2, n = 1001)
s2 <- approxfun(x2, f2,method = nspline_method)
plot(x2, f2, ylim = c(0, 1))
curve(s2(x), add = TRUE, col = 2, n = 1001)
s3 <- function(x){
return( 1 - s1(x) - s2(x) )
}
plot(x3, f3, ylim = c(0, 1))
curve(s3(x), add = TRUE, col = 2, n = 1001)
}
# Plot functions, see if conform sketch
p <- ggplot(data.frame(x=c(0,90)), aes(x))
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s1(x)}, geom="line", colour="red")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)}, geom="line", colour="blue")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)+s3(x)}, geom="line", colour="green")
print (p)
# So, maximize sales of *B* by changing how many features *C* (both at same price)
p <- ggplot(data.frame(x=c(0,90)), aes(x))
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s2(x)}, geom="line", colour="blue")
print (p)
############################## run a test of probs #####################
# set amount of features to 60
feat = 60
choices_sampling <- sample( c("A","B","C"), customers, replace=TRUE, prob=c(s1(feat), s2(feat), s3(feat)))
# plot barchart of prob
choicetable = prop.table(table(choices_sampling))
barplot(choicetable, col = heat.colors(3), ylim=c(0,1))
################################ Do lif! ################################
# now see if LIF can find the optimum C
# which optimizes B, where B is most profitable
yws <- rep(NA, T)
ytm <- .3
length_features <- 90
profit_per_sale <- rep(NA,length_features)
arr.features <- rep(NA, TT)
arr.amount_features_start <- rep(NA,TT)
arr.doesbuy <- rep(NA, TT)
arr.sampleChoices <- rep(NA, TT)
############################## Lif sim loop ################################
for(i in 1:customers){
# Oscillate around features
features <- amount_features_start + A*cos(w*i)
# This is only necessary for non-spline constant
if(features > 90) features <- 90
if(features < 0) features <- 0
sampleChoice <- sample( c("A","B","C"), 1, replace=TRUE, prob=c(s1(features), s2(features), s3(features)))
if (sampleChoice == "A") {
yt <- 30*0.05 # 5% profit on $30 of A
} else if (sampleChoice == "B"){
yt <- 60*0.06 # 6% profit of $60 for B
} else {
yt <- 60*0.04 # 4% profit on $60 for C
}
# rolling average of yt
ytm <- ytm + ((yt - ytm) / window.length)
yws <- push(yws, ytm * cos(w*i))
# reset amount_features_start using lif
if(i > T){
yw <- sum(yws) / T
amount_features_start <- amount_features_start + (gamma / T) * yw
}
# Log me!
profit_per_sale[i] <- yt
arr.sampleChoices[i] <- sampleChoice
arr.features[i] <- features
arr.amount_features_start[i] <- amount_features_start
}
############################## Model profit ###################################
profit_plot <- rep(NA,length_features)
for (j in 1:length_features) profit_plot[j] <- s1(j)*30*0.05 + s2(j)*60*0.06 +s3(j)*60*0.04
############################## Result #########################################
max_profit_model = which(profit_plot==max(profit_plot))
print(paste0("Maximum profit at feature X of profit function: ", max_profit_model))
## [1] "Maximum profit at feature X of profit function: 63"
print(paste0("Maximum profit at feature X found by LiF ", amount_features_start))
## [1] "Maximum profit at feature X found by LiF 59.8611612800204"
############################## do some plotting ##############################
# Do some simple plots
plot(arr.features, type="l",ylim=c(0,90))
plot(profit_plot, type="l")
lines(c(max_profit_model[1],max_profit_model[1]),c(-5,max(profit_plot)), col="red")
plot(slideFunct(profit_per_sale,100,50), type="l")
plot(arr.amount_features_start, type="l", xlab="Time in stream", ylab="features",ylim=c(0,90))
lines(arr.amount_features_start, col="red",ylim=c(0,90))
