We ran the R simulation below to test our hypothesis that LiF should be able to optimize decoy effects (I will add comments further clarifying my train of thought as reflected by the script within the next few days). As the simulation seems to pan out quite well, our next step will be to test our hypothesis in a more realistic setting: online, with real live participants!

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 <- 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 62.7493552775659"

############################## 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))