by Fedor Iskhakov, ANU
Description: Implementation of Rust model in infinite horizon with value function iterations solver
$ EV(x,d) $ = the expected value function
$$ \begin{eqnarray} EV(x,d) &=& \sum_{x' \in X} \log \big( \exp[v(x',0)] + \exp[v(x',1)] \big) \pi(x'|x,d) \\ v(x,d) &=& u(x,d) + \beta EV(x,d) \end{eqnarray} $$Given by the cost function that depends on the choice
$$ u(x,d,\theta_1)=\left \{ \begin{array}{ll} -RC-c(0,\theta_1) & \text{if }d_{t}=\text{replace}=1 \\ -c(x,\theta_1) & \text{if }d_{t}=\text{keep}=0 \end{array} \right. $$If not replacing ($ d=0) $
$$ \Pi(d=0)_{n x n} = \begin{pmatrix} \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ 0 & \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & 0 \\ 0 & 0 &\theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & \cdot & \cdot & 0 & \theta_{20} & \theta_{21} & \theta_{22} & 0 \\ 0 & \cdot & \cdot & \cdot & 0 & \theta_{20} & \theta_{21} & \theta_{22} \\ 0 & \cdot & \cdot & \cdot & \cdot & 0 & \theta_{20} & 1-\theta_{20} \\ 0 & \cdot & \cdot & \cdot & \cdot & \cdot & 0 & 1 \end{pmatrix} $$If replacing ($ d=1) $
$$ \Pi(d=1)_{n x n} = \begin{pmatrix} \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \theta_{20} & \theta_{21} & \theta_{22} & 0 & \cdot & \cdot & \cdot & 0 \\ \end{pmatrix} $$It is always useful to examine the structure of the problem for insights into better code structure!
# THIS IS DEVELOPPED IN THE VIDEO
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
class zurcher():
'''Harold Zurcher bus engine replacement model class, VFI version'''
def __init__(self,
n = 175, # number of state points
RC = 11.7257, # replacement cost
c = 2.45569, # parameter of maintance cost (theta_1)
p = [0.0937,0.4475,0.4459,0.0127], # probabilities of transitions (theta_2)
beta = 0.9999): # discount factor
'''Init for the Zurcher model object'''
assert sum(p)<=1.0, 'Transition probability parameters must sum up to <1'
self.RC, self.c, self.p, self.beta, self.n= RC, c, p, beta, n
@property
def n(self):
'''Attrinute getter for n'''
return self.__n
@n.setter
def n(self, value):
'''Attribute n setter'''
self.__n = value
self.grid = np.arange(self.__n)
self.trpr = self.__transition_probs()
def __repr__(self):
'''String representation of the Zurcher model'''
return 'Rust model at id={}'.format(id(self))
def __transition_probs(self):
'''Computing the transision probability matrix'''
trpr = np.zeros((self.__n,self.__n)) # init
probs = self.p + [1-sum(self.p)] # ensure sum up to 1
for i,p in enumerate(probs):
trpr += np.diag([p]*(self.__n-i),k=i)
trpr[:,-1] = 1.-np.sum(trpr[:,:-1],axis=1)
return trpr
def bellman(self,ev0):
'''Bellman operator for the model'''
x = self.grid # points in the next period state
mcost = -0.001*x*self.c
vx0 = mcost + self.beta * ev0 # vector of v(x,0)
vx1 = -self.RC + mcost[0] + self.beta * ev0[0] # vector of v(x,1)
M = np.maximum(vx0,vx1)
logsum = M + np.log(np.exp(vx0-M) + np.exp(vx1-M))
ev1 = self.trpr @ logsum[:,np.newaxis]
ev1 = ev1.ravel()
# recompute v(x,d) for current period choice probabilities
vx0 = mcost + self.beta * ev1 # vector of v(x,0)
vx1 = -self.RC + mcost[0] + self.beta * ev1[0] # vector of v(x,1)
pr1 = 1 / (np.exp(vx0-vx1)+1) # choice prob to repace
return ev1, pr1
def solve_vfi(self,tol=1e-6,maxiter=100,callback=None):
'''Solves the Rust model using value function iterations
'''
ev0 = np.zeros(self.n) # initial point for VFI
for i in range(maxiter): # main loop
ev1, pr1 = self.bellman(ev0) # update approximation
err = np.amax(np.abs(ev0-ev1))
if callback != None: callback(iter=i,err=err,ev1=ev1,ev0=ev0,pr1=pr1,model=self)
if err<tol:
break # break out if converged
ev0 = ev1 # get ready to the next iteration
else:
raise RuntimeError('Failed to converge in %d iterations'%maxiter)
return ev1, pr1
def solve_show(self,maxiter=1000,tol=1e-6,**kvargs):
'''Illustrate solution'''
fig1, (ax1,ax2) = plt.subplots(1,2,figsize=(14,8))
ax1.grid(b=True, which='both', color='0.65', linestyle='-')
ax2.grid(b=True, which='both', color='0.65', linestyle='-')
ax1.set_xlabel('Mileage grid')
ax2.set_xlabel('Mileage grid')
ax1.set_title('Value function')
ax2.set_title('Probability of replacing the engine')
def callback(**argvars):
mod, ev, pk = argvars['model'],argvars['ev1'],argvars['pr1']
ax1.plot(mod.grid,ev,color='k',alpha=0.25)
ax2.plot(mod.grid,pk,color='k',alpha=0.25)
ev,pk = self.solve_vfi(maxiter=maxiter,tol=tol,callback=callback,**kvargs)
# add solutions
ax1.plot(self.grid,ev,color='r',linewidth=2.5)
ax2.plot(self.grid,pk,color='r',linewidth=2.5)
plt.show()
return ev,pk
model = zurcher(n = 15,p=[0.2,0.2,0.4], beta=0.5)
model = zurcher(beta=0.9)
print(model)
ev,pr = model.solve_show()
print('Solution ev=')
print(ev)
print('Solution pr=')
print(pr)
Rust model at id=140320486659664
Solution ev= [-0.33831227 -0.36286609 -0.38741983 -0.4119735 -0.43652708 -0.46108058 -0.48563399 -0.51018731 -0.53474054 -0.55929368 -0.58384672 -0.60839965 -0.63295248 -0.65750521 -0.68205782 -0.70661032 -0.73116271 -0.75571497 -0.78026711 -0.80481912 -0.829371 -0.85392275 -0.87847435 -0.90302581 -0.92757711 -0.95212827 -0.97667926 -1.00123009 -1.02578075 -1.05033124 -1.07488154 -1.09943165 -1.12398157 -1.14853129 -1.1730808 -1.19763009 -1.22217915 -1.24672799 -1.27127658 -1.29582492 -1.32037301 -1.34492082 -1.36946835 -1.39401558 -1.41856251 -1.44310913 -1.46765541 -1.49220135 -1.51674693 -1.54129213 -1.56583694 -1.59038134 -1.61492531 -1.63946883 -1.66401189 -1.68855444 -1.71309649 -1.73763799 -1.76217892 -1.78671925 -1.81125896 -1.835798 -1.86033635 -1.88487396 -1.9094108 -1.93394683 -1.958482 -1.98301625 -2.00754955 -2.03208184 -2.05661305 -2.08114312 -2.10567199 -2.13019959 -2.15472583 -2.17925063 -2.20377391 -2.22829557 -2.25281551 -2.27733362 -2.30184978 -2.32636387 -2.35087575 -2.37538527 -2.39989229 -2.42439663 -2.44889812 -2.47339655 -2.49789174 -2.52238344 -2.54687143 -2.57135545 -2.59583522 -2.62031044 -2.64478079 -2.66924593 -2.6937055 -2.71815908 -2.74260625 -2.76704656 -2.7914795 -2.81590455 -2.84032112 -2.8647286 -2.88912631 -2.91351355 -2.93788954 -2.96225345 -2.98660438 -3.01094136 -3.03526335 -3.05956923 -3.0838578 -3.10812774 -3.13237767 -3.15660606 -3.1808113 -3.20499161 -3.22914512 -3.25326977 -3.27736338 -3.30142358 -3.3254478 -3.3494333 -3.37337711 -3.39727603 -3.4211266 -3.44492512 -3.46866756 -3.4923496 -3.51596659 -3.53951349 -3.56298487 -3.58637491 -3.60967728 -3.63288518 -3.6559913 -3.67898771 -3.70186589 -3.72461663 -3.74723003 -3.76969536 -3.7920011 -3.81413479 -3.836083 -3.85783125 -3.87936391 -3.90066412 -3.92171371 -3.94249306 -3.962981 -3.98315473 -4.00298961 -4.02245908 -4.04153448 -4.0601849 -4.07837698 -4.09607474 -4.11323936 -4.12982895 -4.14579834 -4.16109875 -4.17567765 -4.18947815 -4.20243932 -4.21449448 -4.22557372 -4.23559719 -4.24448843 -4.25214374 -4.25849757 -4.2633781 -4.26682051 -4.26838168 -4.26863285] Solution pr= [8.08331835e-06 8.28425234e-06 8.49018048e-06 8.70122688e-06 8.91751870e-06 9.13918627e-06 9.36636317e-06 9.59918628e-06 9.83779588e-06 1.00823357e-05 1.03329532e-05 1.05897993e-05 1.08530287e-05 1.11228000e-05 1.13992757e-05 1.16826225e-05 1.19730109e-05 1.22706159e-05 1.25756168e-05 1.28881972e-05 1.32085456e-05 1.35368547e-05 1.38733223e-05 1.42181511e-05 1.45715487e-05 1.49337280e-05 1.53049069e-05 1.56853091e-05 1.60751635e-05 1.64747049e-05 1.68841737e-05 1.73038166e-05 1.77338861e-05 1.81746410e-05 1.86263467e-05 1.90892750e-05 1.95637044e-05 2.00499205e-05 2.05482158e-05 2.10588901e-05 2.15822505e-05 2.21186120e-05 2.26682970e-05 2.32316363e-05 2.38089684e-05 2.44006405e-05 2.50070084e-05 2.56284365e-05 2.62652982e-05 2.69179763e-05 2.75868629e-05 2.82723598e-05 2.89748787e-05 2.96948415e-05 3.04326803e-05 3.11888381e-05 3.19637687e-05 3.27579368e-05 3.35718189e-05 3.44059030e-05 3.52606890e-05 3.61366892e-05 3.70344283e-05 3.79544440e-05 3.88972869e-05 3.98635212e-05 4.08537248e-05 4.18684897e-05 4.29084221e-05 4.39741432e-05 4.50662890e-05 4.61855109e-05 4.73324761e-05 4.85078678e-05 4.97123856e-05 5.09467459e-05 5.22116822e-05 5.35079454e-05 5.48363043e-05 5.61975458e-05 5.75924756e-05 5.90219180e-05 6.04867167e-05 6.19877351e-05 6.35258564e-05 6.51019844e-05 6.67170434e-05 6.83719787e-05 7.00677571e-05 7.18053669e-05 7.35858183e-05 7.54101440e-05 7.72793990e-05 7.91946610e-05 8.11570309e-05 8.31676324e-05 8.52276126e-05 8.73381422e-05 8.95004150e-05 9.17156485e-05 9.39850834e-05 9.63099838e-05 9.86916368e-05 1.01131352e-04 1.03630462e-04 1.06190322e-04 1.08812305e-04 1.11497810e-04 1.14248251e-04 1.17065063e-04 1.19949698e-04 1.22903624e-04 1.25928323e-04 1.29025289e-04 1.32196027e-04 1.35442050e-04 1.38764875e-04 1.42166019e-04 1.45647000e-04 1.49209328e-04 1.52854501e-04 1.56584004e-04 1.60399296e-04 1.64301811e-04 1.68292943e-04 1.72374046e-04 1.76546419e-04 1.80811296e-04 1.85169839e-04 1.89623121e-04 1.94172115e-04 1.98817677e-04 2.03560530e-04 2.08401244e-04 2.13340219e-04 2.18377656e-04 2.23513539e-04 2.28747600e-04 2.34079296e-04 2.39507771e-04 2.45031820e-04 2.50649852e-04 2.56359843e-04 2.62159292e-04 2.68045167e-04 2.74013847e-04 2.80061067e-04 2.86181846e-04 2.92370421e-04 2.98620165e-04 3.04923508e-04 3.11271849e-04 3.17655461e-04 3.24063389e-04 3.30483347e-04 3.36901605e-04 3.43302875e-04 3.49670185e-04 3.55984757e-04 3.62225879e-04 3.68370781e-04 3.74394488e-04 3.80269739e-04 3.85966772e-04 3.91453422e-04 3.96694563e-04 4.01653058e-04 4.06287423e-04 4.10556841e-04 4.14410073e-04 4.17809980e-04 4.20679823e-04 4.23021728e-04 4.24657367e-04 4.25797233e-04]
# THIS IS AN EARLIER VERSION OF THE CODE
# (there is bug here..)
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
class zurcher():
'''Harold Zurcher bus engine replacement model class, VFI version'''
def __init__(self,
n = 175, # number of state points
RC = 11.7257, # replacement cost
c = 2.45569, # parameter of maintance cost (theta_1)
p = [0.0937,0.4475,0.4459,0.0127], # probabilities of transisions (theta_2)
beta = 0.9999): # discount factor
'''Initializator for the zurcher class'''
assert sum(p)<=1.0,'sum of transision probability parameters must not exceed unity'
self.p = p # parameters for transision probabilities
self.n = n # set number of grid points on the state space
self.RC = RC # replacement cost
self.c = c # cost function parameter
self.beta = beta # discount factor
@property
def n(self):
'''Dimension getter'''
return self.__n # internal dimension variable
@n.setter
def n(self,n):
'''Dimension setter, updaing the grid and transision probabilities'''
self.__n = n
# create grid for the set dimension
self.grid = np.arange(n) # 0,..,n-1 index of grid points
# create transition prob for the set dimension
p = self.p # "copy" the list of parameters
p.append(1.0-sum(p)) # add the last element to ensure 1.0 in sum
self.P1,self.P2 = self.transition_probability(np.array(p)) # compute transision probabilities
def __repr__(self):
'''String representation of the Zurcher model object'''
# id() is unique identifier for the variable (reference), convert to hex
return 'Zurcher bus engine replacement model with id=%s' % hex(id(self))
def transition_probability(self,p):
'''Compute transition probability matrixes conditional on choice'''
# conditional on d=0, no replacement
P1 = np.full((self.n,self.n),0.0)
for i in range(self.n):
if i <= self.n-p.size:
# lines where p vector fits entirely
P1[i][i:i+p.size]=p
else:
P1[i][i:] = p[:self.n-p.size-i]
P1[i][-1] = 1.0-P1[i][:-1].sum()
# conditional on d=1, replacement
P2 = np.full((self.n,self.n),0.0)
for i in range(self.n):
P2[i][:p.size]=p
return P1,P2
def bellman(self,ev0,output=0):
''' Bellman operator for the model
Input: current approximation of the EV as column vector
output = type of output requested
Output: new approximation of EV
d=0 choice probability (if output>0)
Frechet derivative of Bellman operator (if output>1)
'''
# EV0 is current approximation of EV on the fixed grid
# For d=0 it holds values for all mileages
# For d=1 (replacement) we use the first value EV0[0]
# So, instead of EV(x,d) for d=0,1, we can use only one vector!
assert np.all(ev0.shape==(self.n,1)),'Expecting EV as column vector'
x = self.grid.reshape((self.n,1)) # states (in the next period), column vector
c = 0.001*self.c*x # maintenance cost in all states
v0 = -c + self.beta*ev0 # value of not replacing
v1 = -c[0] -self.RC + self.beta*ev0[0] # value of replacing
# recenter the values for numerical stability of logsum !!!!!!!!!!!!!!!!
maxv = np.maximum(v0,v1)
logsum = maxv + np.log(np.exp(v0-maxv) + np.exp(v1-maxv))
ev1 = self.P1 @ logsum # matrix multiplication, result as column vector
if output == 0:
return ev1
# keep (no replacement) choice probability
pk = 1/(1+np.exp(v1-v0))
if output == 1:
return ev1,pk
def solve_vfi(self, maxiter=1000, tol=1e-6, callback=None):
'''Solves the model using successive approximations (VFI)'''
ev0 = np.full((self.n,1),0) # initial guess of EV
for iter in range(maxiter):
ev1,pk = self.bellman(ev0,output=1)
stp = np.max(abs(ev1-ev0))
if callback:
if iter==0: stp0=1.0
callback(iter,self,ev1,pk,stp,stp/stp0) # callback for making plots
if stp < tol:
break
ev0=ev1
stp0=stp
else: # when i went up to maxiter
print('No convergence: maximum number of iterations achieved!')
return ev1,pk
def solve_show(self,solver='vfi',maxiter=1000,tol=1e-6,**kvargs):
'''Illustrate solution'''
fig1, (ax1,ax2) = plt.subplots(1,2,figsize=(14,8))
ax1.grid(b=True, which='both', color='0.65', linestyle='-')
ax2.grid(b=True, which='both', color='0.65', linestyle='-')
ax1.set_xlabel('Mileage grid')
ax2.set_xlabel('Mileage grid')
ax1.set_title('Value function')
ax2.set_title('Probability of keeping the engine')
def callback(iter,mod,ev,pk,stp,dstp):
if iter==0:
print('%4s %16s %16s'%('iter','err','err(i)/err(i-1)'))
print('-'*40)
print('%4d %16.12f %16.12f'%(iter,stp,dstp))
ax1.plot(mod.grid,ev,color='k',alpha=0.25)
ax2.plot(mod.grid,pk,color='k',alpha=0.25)
if solver=='vfi':
ev,pk = self.solve_vfi(maxiter=maxiter,tol=tol,callback=callback,**kvargs)
elif solver=='nk':
ev,pk = self.solve_nk(maxiter=maxiter,tol=tol,callback=callback,**kvargs)
elif solver=='poly':
ev,pk = self.solve_poly(maxiter=maxiter,tol=tol,callback=callback,**kvargs)
else:
print('Unknown solver')
return None
# add solutions
ax1.plot(self.grid,ev,color='r',linewidth=2.5)
ax2.plot(self.grid,pk,color='r',linewidth=2.5)
plt.show()
return ev,pk
# investigate how parts of the code work:
model = zurcher(RC=.5,n=12,p=[0.65,0.2,0.1]) # model instance
print('Model grid:\n',model.grid)
print(model) # string representation
print('Transition probabilities conditional on not replacing:\n',model.P1)
print('Transition probabilities conditional on replacing:\n',model.P2)
ev,pk = model.bellman(np.full((model.n,1),0),output=1)
print('Bellman one run:\n',ev)
print('Probability of keeping:\n',pk)
Model grid: [ 0 1 2 3 4 5 6 7 8 9 10 11] Zurcher bus engine replacement model with id=0x7f9eb83bda90 Transition probabilities conditional on not replacing: [[0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0.65 0.2 0.1 0.05 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0.65 0.2 0.1 0.05 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0.65 0.2 0.1 0.05 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0.65 0.2 0.1 0.05 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0.65 0.2 0.1 0.05] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.65 0.2 0.15] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.65 0.35] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. ]] Transition probabilities conditional on replacing: [[0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ] [0.65 0.2 0.1 0.05 0. 0. 0. 0. 0. 0. 0. 0. ]] Bellman one run: [[0.47323702] [0.47170994] [0.47018428] [0.46866004] [0.46713722] [0.46561582] [0.46409584] [0.46257729] [0.46106015] [0.45962006] [0.45833254] [0.45734867]] Probability of keeping: [[0.62245933] [0.62188206] [0.62130445] [0.62072649] [0.62014818] [0.61956954] [0.61899056] [0.61841124] [0.61783158] [0.61725158] [0.61667125] [0.61609059]]
# solve Harold Zurcher model for different parameters
m = zurcher(beta=0.9,RC=5.0)
ev,pk = m.solve_show(maxiter=500,tol=1e-4)
iter err err(i)/err(i-1) ---------------------------------------- 0 0.417013120838 0.417013120838 1 0.370621231082 0.888751966213 2 0.327694171617 0.884175390224 3 0.287961455490 0.878750616984 4 0.251303448761 0.872698216966 5 0.217734713293 0.866421509000 6 0.187291798495 0.860183457487 7 0.160016330707 0.854369128773 8 0.135901357738 0.849296800756 9 0.114862847118 0.845192785634 10 0.096741182600 0.842232149276 11 0.081326735019 0.840663023065 12 0.068378239873 0.840784274152 13 0.057792280046 0.845185254165 14 0.049175736163 0.850904932690 15 0.042089382417 0.855897353068 16 0.036210267283 0.860318332170 17 0.031292801483 0.864196920689 18 0.027151337854 0.867654430650 19 0.023641549056 0.870732380965 20 0.020650667367 0.873490451829 21 0.018089047165 0.875954604423 22 0.015885781951 0.878198934777 23 0.013982833446 0.880210586351 24 0.012332999988 0.882010075810 25 0.010898109646 0.883654395295 26 0.009646543045 0.885157459273 27 0.008551907007 0.886525563278 28 0.007592036396 0.887759465807 29 0.006748259343 0.888860246554 30 0.006005359648 0.889912397087 31 0.005349690415 0.890819322886 32 0.004770277908 0.891692329344 33 0.004257352714 0.892474777412 34 0.003802587451 0.893181210732 35 0.003398868436 0.893830445683 36 0.003040032902 0.894425000077 37 0.002720730799 0.894967550499 38 0.002436307958 0.895460865964 39 0.002182707191 0.895907754115 40 0.001956454399 0.896343039906 41 0.001754395785 0.896722042706 42 0.001573785710 0.897052833441 43 0.001412253499 0.897360733435 44 0.001267744249 0.897674709159 45 0.001138319591 0.897909489115 46 0.001022391934 0.898158954578 47 0.000918488806 0.898372508819 48 0.000825313079 0.898555402383 49 0.000741749231 0.898748910814 50 0.000666747945 0.898885927178 51 0.000599439885 0.899050216880 52 0.000538991569 0.899158669047 53 0.000484708110 0.899286998748 54 0.000435937364 0.899381204954 55 0.000392112077 0.899468843198 56 0.000352728537 0.899560502761 57 0.000317316129 0.899604358475 58 0.000285490244 0.899702906408 59 0.000256863810 0.899728853221 60 0.000231122587 0.899786493992 61 0.000207970483 0.899827604103 62 0.000187139721 0.899837888320 63 0.000168407479 0.899902375562 64 0.000151550220 0.899901957177 65 0.000136383280 0.899921362375 66 0.000122738027 0.899949221210 67 0.000110456674 0.899938484602 68 0.000099407298 0.899966423032