1719.sagews

#-*- coding: utf-8 -*- 
# Course: CS250 - Computational Methods 
# Date: 03/04/10 
# Purpose: Recreating Graph 3 from Congdon 
# CODING DONE IN PYTHON! 


from scipy import *
from pylab import *
from math import *
from numpy import *
import matplotlib.pyplot as plt


# upload file chse.txt
file = open ("chse.txt", "r")
Chse = file.readlines() 

for i in xrange(len(Chse)):  
    Chse[i]=Chse[i].split() 
    
for i in range(len(Chse)): 
    for j in range(len(Chse[i])): 
        Chse[i][j]=float(Chse[i][j]) 
        
Age = range(len(Chse))
mx = [Chse[F][3] for F in range (len(Chse))] 
lx = [Chse[p][2] for p in range (len(Chse))]
sx = [Chse[p][1] for p in range (len(Chse))]

Age = array(Age, float)
mx = array(mx, float)
lx = array(lx, float)
sx = array(sx, float)

# Lx list dependent on sx 
a = 1.
Lx = [a]
for k in range (len(lx)):
    a = a * sx[k]
    Lx.append(a.tolist())

# creating sx arrays
Px = array([0.01, 0.2, 0.4, 0.6, 0.8, 1.])
sx0 = array(sx,float)
sx2 = array(sx,float)
sx4 = array(sx,float)
sx6 = array(sx,float)
sx8 = array(sx,float)
sx10 = array(sx,float)

sx0[1:12:1] = Px[0]
sx2[1:12:1] = Px[1]
sx4[1:12:1] = Px[2]
sx6[1:12:1] = Px[3]
sx8[1:12:1] = Px[4]
sx10[1:12:1] = Px[5]

# creating Lx arrays 
a0 = 1.
Lx0 = [a0]
for l in range (len(lx)-1):
    a0 = a0 * sx0[l]
    Lx0.append(a0.tolist())
Lx0 = array(Lx0, float)

a2 = 1.
Lx2 = [a2]
for k in range (len(lx)-1):
    a2 = a2 * sx2[k]
    Lx2.append(a2.tolist())
Lx2 = array(Lx2, float)

a4 = 1.
Lx4 = [a4]
for k in range (len(lx)-1):
    a4 = a4 * sx4[k]
    Lx4.append(a4.tolist())
Lx4 = array(Lx4, float)

a6 = 1.
Lx6 = [a6]
for k in range (len(lx)-1):
    a6 = a6 * sx6[k]
    Lx6.append(a6.tolist())
Lx6 = array(Lx6, float)

a8 = 1.
Lx8 = [a8]
for k in range (len(lx)-1):
    a8 = a8 * sx8[k]
    Lx8.append(a8.tolist())
Lx8 = array(Lx8, float)

a10 = 1.
Lx10 = [a10]
for k in range (len(lx)-1):
    a10 = a10 * sx10[k]
    Lx10.append(a10.tolist())
Lx10 = array(Lx10, float)

# r values

def calc_error(lxmx, r, Age):
    return 1. - sum(lxmx * e**(-r*Age))

lxmx = lx * mx
lxmxx = lx * mx * Age

r0 = sum(lxmx)
G = sum(lxmxx) / sum(lxmx)
r_est = log(r0)/G 

error = calc_error(lxmx, r_est, Age)
tolerance = 0.0000000000001

if error >= 0:
    r_low = r_est - 1.0
    r_high = r_est 
else:
    r_low = r_est
    r_high = r_est + 1.0
    
while abs(error) >= tolerance:
    r = (r_low + r_high) / 2.
    error = calc_error(lxmx, r, Age)
    if error >= 0:
        r_high = r
    else:
        r_low = r

# r0

lxmx0 = Lx0 * mx
lxmxx0 = Lx0 * mx * Age

r0_0 = sum(lxmx)
G0 = sum(lxmxx) / sum(lxmx)
r_est0 = log(r0_0)/G0 

error0 = calc_error(lxmx0, r_est0, Age)
tolerance = 0.0000000000001

if error >= 0:
    r_low = r_est0 - 1.0
    r_high = r_est0 
else:
    r_low = r_est0
    r_high = r_est0 + 1.0
    
while abs(error0) >= tolerance:
    r0 = (r_low + r_high) / 2.
    error0 = calc_error(lxmx0, r0, Age)
    if error0 >= 0:
        r_high = r0
    else:
        r_low = r0

#r2
lxmx2 = Lx2 * mx
lxmxx2 = Lx2 * mx * Age

r0_2 = sum(lxmx2)
G2 = sum(lxmxx2) / sum(lxmx2)
r_est2 = log(r0_2)/G2 

error2 = calc_error(lxmx2, r_est2, Age)
tolerance = 0.0000000000001

if error2 >= 0:
    r_low = r_est2 - 1.0
    r_high = r_est2 
else:
    r_low = r_est2
    r_high = r_est2 + 1.0
    
while abs(error2) >= tolerance:
    r2 = (r_low + r_high) / 2.
    error2 = calc_error(lxmx2, r2, Age)
    if error2 >= 0:
        r_high = r2
    else:
        r_low = r2

# r4
lxmx4 = Lx4 * mx
lxmxx4 = Lx4 * mx * Age

r0_4 = sum(lxmx4)
G4 = sum(lxmxx4) / sum(lxmx4)
r_est4 = log(r0_4)/G4 

error4 = calc_error(lxmx4, r_est4, Age)
tolerance = 0.0000000000001

if error4 >= 0:
    r_low = r_est4 - 1.0
    r_high = r_est4 
else:
    r_low = r_est4
    r_high = r_est4 + 1.0
    
while abs(error4) >= tolerance:
    r4 = (r_low + r_high) / 2.
    error4 = calc_error(lxmx4, r4, Age)
    if error4 >= 0:
        r_high = r4
    else:
        r_low = r4

#r6
lxmx6 = Lx6 * mx
lxmxx6 = Lx6 * mx * Age

r0_6 = sum(lxmx6)
G6 = sum(lxmxx6) / sum(lxmx6)
r_est6 = log(r0_6)/G6 

error6 = calc_error(lxmx6, r_est6, Age)
tolerance = 0.0000000000001

if error6 >= 0:
    r_low = r_est6 - 1.0
    r_high = r_est6 
else:
    r_low = r_est6
    r_high = r_est6 + 1.0
    
while abs(error6) >= tolerance:
    r6 = (r_low + r_high) / 2.
    error6 = calc_error(lxmx6, r6, Age)
    if error6 >= 0:
        r_high = r6
    else:
        r_low = r6

#r8
lxmx8 = Lx8 * mx
lxmxx8 = Lx8 * mx * Age

r0_8 = sum(lxmx8)
G8 = sum(lxmxx8) / sum(lxmx8)
r_est8 = log(r0_8)/G8

error8 = calc_error(lxmx8, r_est8, Age)
tolerance = 0.0000000000001

if error8 >= 0:
    r_low = r_est8 - 1.0
    r_high = r_est8 
else:
    r_low = r_est8
    r_high = r_est8 + 1.0
    
while abs(error8) >= tolerance:
    r8 = (r_low + r_high) / 2.
    error8 = calc_error(lxmx8, r8, Age)
    if error8 >= 0:
        r_high = r8
    else:
        r_low = r8

# r10
lxmx10 = Lx10 * mx
lxmxx10 = Lx10 * mx * Age

r0_10 = sum(lxmx10)
G10 = sum(lxmxx10) / sum(lxmx10)
r_est10 = log(r0_10)/G10 

error10 = calc_error(lxmx10, r_est10, Age)
tolerance = 0.0000000000001

if error10 >= 0:
    r_low = r_est10 - 1.0
    r_high = r_est10 
else:
    r_low = r_est10
    r_high = r_est10 + 1.0
    
while abs(error10) >= tolerance:
    r10 = (r_low + r_high) / 2.
    error10 = calc_error(lxmx10, r10, Age)
    if error10 >= 0:
        r_high = r10
    else:
        r_low = r10

#print r0, r2, r4, r6, r8, r10

x_values = [0.01, 0.2, 0.4, 0.6, 0.8, 1.]

r_values = [r0, r2, r4, r6, r8, r10]


matplotlib.pyplot.axhline(linewidth=1, color='black')
plt.plot( x_values, r_values)
plt.show()