In [1]:
#library
import numpy as np
import math 
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import pearsonr as pcoor

#Data
Amount = 3000
Delays = [7, 30, 180, 365, 1085] 
IndifferencePoints = [2800, 1850, 800, 850, 550]

#Code
z = np.full(len(Delays), Amount)
x = np.array(Delays)
y = np.array(IndifferencePoints)

def Rely(x,y):
    return x / max(y) 

def Eq0(x, k): 
    return z*np.exp(k*-x)
kEq0, _ = curve_fit(Eq0, x, y, maxfev = 40000) 
yEq0 = Eq0(x,kEq0)
rEq0, pEq0 = pcoor(y, yEq0)

def Eq1(x, k): 
    return z/(1+(x*k))
kEq1, _ = curve_fit(Eq1, x, y, maxfev = 40000) 
yEq1 = Eq1(x,kEq1)
rEq1, pEq1 = pcoor(y, yEq1)

def Eq2(x, k, s): 
    return z/((1+(x*k))*np.exp(s))
kEq2, _ = curve_fit(Eq2, x, y, maxfev = 40000) 
yEq2 = Eq2(x,kEq2[0],kEq2[1])
rEq2, pEq2 = pcoor(y, yEq2)

def Eq3(x, k, s): 
    return z/(1+(x*(k*np.exp(s))))
kEq3, _ = curve_fit(Eq3, x, y, maxfev = 40000) 
yEq3 = Eq3(x,kEq3[0],kEq3[1])
rEq3, pEq3 = pcoor(y, yEq3)

def AuC(x,y):
    temp = 0
    tempr = 0
    for i in range(1,len(x)):
        temp += (min(y[i],y[i-1])*abs(x[i]-x[i-1]))+((abs(max(y[i],y[i-1])-min(y[i],y[i-1]))*abs(x[i]-x[i-1]))/2)
        tempr += (min(Rely(y,z)[i],Rely(y,z)[i-1])*abs(x[i]-x[i-1]))+((abs(max(Rely(y,z)[i],Rely(y,z)[i-1])-min(Rely(y,z)[i],Rely(y,z)[i-1]))*abs(x[i]-x[i-1]))/2)
    return math.sqrt(temp), math.sqrt(tempr) 
    
plt.figure(1)    
plt.plot(x, Rely(yEq0,z), 'k--', label='Exponential') 
plt.plot(x, Rely(yEq1,z), 'r--', label='Mazur') 
plt.plot(x, Rely(yEq2,z), 'g--', label='Myerson and Green') 
#plt.plot(x, Rely(yEq3,z), 'b--', label='Rachlin') 
plt.plot(x, Rely(y,z), 'ko') 
plt.xlabel('Delay')
plt.ylabel('Indifference Points')
plt.legend()
plt.axis([min(x)-(max(x)*.1), max(x)+(max(x)*.1), 0, 1])

print('Exponential      :','R²=' + str(round(rEq0,3)), 'k='+str(round(kEq0[0],3)))
print('Mazur            :','R²=' + str(round(rEq1,3)), 'k='+str(round(kEq1[0],3)))
print('Myerson and Green:','R²=' + str(round(rEq2,3)), 'k='+str(round(kEq2[0],3)), 's='+str(round(kEq2[1],3)))
print('Rachlin          :','R²=' + str(round(rEq3,3)), 'k='+str(round(kEq3[0],3)), 's='+str(round(kEq3[1],3)))
print('')
print('Area Under Curve:','AuC²=' + str(round(AuC(x,y)[0],3)))
print('Relative AuC    :','AuC²=' + str(round(AuC(x,y)[1],3)))
Exponential      : R²=0.951 k=0.007
Mazur            : R²=0.975 k=0.013
Myerson and Green: R²=0.967 k=0.011 s=0.06
Rachlin          : R²=0.975 k=0.054 s=-1.38

Area Under Curve: AuC²=953.336
Relative AuC    : AuC²=17.405