LiF_Decoy_Local_Test.py

# -*- coding: utf-8 -*-

#==============================================================================

# Reset IPython - may delete the following when not making use of Spyder IDE

from IPython import get_ipython
get_ipython().magic('reset -sf')

# import libraries

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
from scipy import misc

#==============================================================================

def matrixpush(m, row):
    if not np.all(np.isfinite(values[:,0])):
        i = np.count_nonzero(np.logical_not(np.isnan(values[:,0])))
        m[i,] = row
    else:
        m = np.vstack([m,row])
        m = m[1:,]
    return(m)

#==============================================================================


x_value_min = 0                   # Domain, min X 
x_value_max = 1                   # Domain, max X, is of influence on 
                                  # Start value and *Amplitude* !

x0 = 0.3 * x_value_max            # Start value

rounded_to = 1                    # X shown to subject is rounded
                                  # to the given number of decimals. 
                                  # Rounded_to has to round to at least 
                                  # an order of magnitude divided by 2 
                                  # smaller values than x_value_max.
                                  # If not, the displayed x
                                  # would end up being 0.0 all of the time.

y_line = 850                      # y of decoy line

p_return = 0.95                   # Probability that agent makes a choice

stream = 3000                     # Length of stream
inttime = 150                     # Integration time
amplitude = 0.07 * x_value_max    # Amlitude
learnrate = .06                   # Learnrate
omega = 1.0                       # Omega

#==============================================================================

# initialize some variables

values = np.zeros((inttime,3))   # prefill value matrix
values.fill(np.nan)              # with np.nan

track_x0 = []                    # initialize x0 log
track_x  = []                    # initialize x log

x = 0.0                          # set x to zero
t = 0.0                          # set t to zero
y = 0.0                          # set y to zero

#==============================================================================
#
# In the current simulation we make use of a Decoy decision model from:
#
# Trueblood, Jennifer S., Scott D. Brown, and Andrew Heathcote. 
# "The multiattribute linear ballistic accumulator model of context effects in 
# multialternative choice." 
# Psychological review 121.2 (2014): 179.
#
# The model of Trueblood et al is itself based on data from previous
# experiments. We will use it here to set the parameters of our simulation.
#
#==============================================================================

# import matrix graph from article

decoy_img = misc.imread("Decoy_Matrix.bmp")

plt.imshow(decoy_img)
plt.axhline(y=y_line,linewidth=1, color='g',linestyle='--')
plt.text(650-20, 650+20, "A", fontsize=12, color='w')
plt.text(decoy_img.shape[0]-650-20, 
         decoy_img.shape[0]-650-20, "B", fontsize=12, color='black')
         
plt.show()

# create brightness matrix 

decoy_img_brightness = decoy_img.sum(axis=2) / decoy_img.shape[2]

# set decoy accroding to brightness along the chosen line
# approximation, brightness seems not to be 1:1 to p, 
# but good enough for our purpose

decoy_y = decoy_img_brightness[y_line] / 255

# create xvalues, from x min to x max

decoy_x = np.arange(x_value_min, decoy_y.size)
decoy_x = (decoy_x / decoy_y.size) * x_value_max

# plot the decoy values as based on brightness data

plt.plot(decoy_x, decoy_y,linewidth=1,color="blue",linestyle='-')


# create a spline approximation of the decoy values 
# based on brightness data

spline = UnivariateSpline (decoy_x, decoy_y-0.04, k=4, s=1)
spline_roots = spline.derivative().roots()

# relevant local max - hardcoded, 
# might be neat to use second deriv pos check, ie
# spline.derivatives(spline.derivative().roots()[x])[1] == postive

spline_local_max = spline_roots[1]
plt.axvline(x=spline_local_max,linewidth=1, color='r')

# plot the spline in the same graph

plt.plot(decoy_x,spline(decoy_x), color="green",linewidth=1)

plt.show()

#==============================================================================
#
# We now get to the main loop of our simulation. It will search for the 
# optimal decoy value, making use of the Lock in Feedback algorithm
# as described in:
# 
# Kaptein, Maurits, and Davide Ianuzzi. 
# "Lock in Feedback in Sequential Experiment." 
# arXiv preprint arXiv:1502.00598 (2015).
# 
#==============================================================================

for i in range(0,stream):

    # set the time

    t = i
    
    # do LiF
    
    x =x0 + amplitude*np.cos(omega * t)
    if np.all(np.isfinite(values[:,0])):
        x0 = np.mean(values[:,1])
        x0 = x0 + learnrate * sum( values[:,2] )
        values.fill(np.nan)


    # just a failsafe, makes sure value not out of bounds, mostly superfluous

    if(x > x_value_max): x = x_value_max
    if(x < x_value_min): x = x_value_min
 
    # in our experiment, displayed x will often be rounded to integers
 
    display_x = np.round(x,rounded_to)

    # randomly pick A's or B's - p changes according to position of decoy C 
    # (funct. np.random.binomial would have been other option for same result)

    sampleChoice = np.random.choice(["B","A"], 
                                    size=(1,), 
                                    p=[spline(display_x), 1-spline(display_x)])  
        
    # did agent/subject choose A or B 
        
    if sampleChoice == 'B':
       yt = 1.0                                  # chooses B! Hurray!
    else:
       yt = 0.0                                  # did not choose B
       
    # does not always make choice   

    if np.random.binomial(1, p_return, 1)==1:    

       y = amplitude*np.cos(omega * t)*yt        # but if did make a choice, 
       row_to_add = np.array([t,x,y])            # enter y into second part LiF
       values = matrixpush(values, row_to_add)   # v1 (set line comment for v2)
       
    # track x and x0 values   

    track_x  = np.append(track_x, x)             # log x
    track_x0 = np.append(track_x0, x0)           # log x0


#==============================================================================

# plot log of x0

plt.plot(track_x0, color='black',linewidth=1)
plt.axhline(y=spline_local_max,linewidth=1, color='r',linestyle='--')
plt.ylim([0,x_value_max/2])
plt.show()

# print last x0 as estimated by LiF

print("Found max: " + str(x0))

# compare to our models max

print("Spline max: "+ str(spline_local_max))


#==============================================================================