Problem 1

Part 3

% matplotlib inline
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
from exponential import *

expodata = np.loadtxt('./expodata.txt',skiprows=1)#first row is the label

#linearize the data
#X(t) = a * t + log(Xo)
expo_log = np.zeros((expodata.shape[0],expodata.shape[1]))
expo_log[:,0] = expodata[:,0]
expo_log[:,1] = np.log(expodata[:,1])
[slope,intercept,r_regres,p_regres,sterr] = stats.linregress(expo_log)

expofit = exponential(expodata[:,0],np.exp(intercept),slope)

plt.plot(expodata[:,0],expodata[:,1],'ro')
plt.plot(expodata[:,0],expofit)
plt.ylabel('X [ind]')
plt.xlabel('t [yrs]')
plt.show()

Part 4

% matplotlib inline
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
from exponential import *
from scipy import optimize

expodata = np.loadtxt('./expodata.txt',skiprows=1)#first row is the label

b,pcov = optimize.curve_fit(exponential,expodata[:,0],expodata[:,1])

expofit = exponential(expodata[:,0],b[0],b[1])

#doubling time is independent of time
doubling_time = np.log(2.0)/b[1]
print ('The time for doubling is estimated as: %1.2e' %(doubling_time)+' yrs') 
#reasons why not: 1: there could be limited carrying capacity, wrong model, extrapolation
#estimated time fro final population to double.
print ('Time to double of final: %1.4e' %(np.log(expodata[10,1]*2./b[0])/b[1]) +' yrs') 

plt.plot(expodata[:,0],expodata[:,1],'ro')
plt.plot(expodata[:,0],expofit)
plt.ylabel('X [ind]')
plt.xlabel('t [yrs]')
plt.show()

Problem 2

% matplotlib inline
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
from exponential import *

expodata = np.loadtxt('./expodata2.txt',skiprows=1)#first row is the label

#Linear Regression
expo_log = np.zeros((expodata.shape[0],expodata.shape[1]))
expo_log[:,0] = expodata[:,0]
expo_log[:,1] = np.log(expodata[:,1])
[slope,intercept,r_regres,p_regres,sterr] = stats.linregress(expo_log)
expofit_l = exponential(expodata[:,0],np.exp(intercept),slope)

#NonLinear Regression
def eqn(x, XO,a):
        return XO * np.exp(a*x)
b_nl,pcov = optimize.curve_fit(exponential,expodata[:,0],expodata[:,1],p0=[1.,-.01])
expofit_nl = exponential(expodata[:,0],b_nl[0],b_nl[1])

#doubling time is independent of time
half_life = -np.log(2.0)/slope
print ('Linear Regression: The half-life is estimated as: %1.4e' %(half_life)+' days') 
print ('Time to half of initial: %1.4e' %(np.log(expodata[0,1]/2./np.exp(intercept))/slope) +' days')
half_life = -np.log(2.0)/b_nl[1]
print ('NonLinear Regression: The half-life is estimated as: %1.4e' %(half_life)+' days')
print ('Time to half of initial: %1.4e' %(np.log(expodata[0,1]/2./b_nl[0])/b_nl[1]) +' days')

fig = plt.figure(1)
fig.subplots_adjust(left=-.5)
plt.subplot(1,2,1)
plt.title('Linear Regression')
plt.plot(expodata[:,0],expodata[:,1],'ro')
plt.plot(expodata[:,0],expofit_l)
plt.ylabel('X [Bq]')
plt.xlabel('t [days]')
plt.subplot(1,2,2)
plt.title('NonLinear Regression')
plt.plot(expodata[:,0],expodata[:,1],'ro')
plt.plot(expodata[:,0],expofit_nl)
plt.ylabel('X [Bq]')
plt.xlabel('t [days]')
plt.show()