# -*- coding: utf-8 -*-
"""
Created on Tue Apr  7 18:40:12 2020

@author: caleb
"""

## Section 1 - Import necessary libraries
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm


## Section 2 - Download and read data file
url = 'https://docs.google.com/spreadsheets/d/e/2PACX-1vRQ6XSqSc79HlLHsOkFHlODwJ7CgHZ8Xp-TegHZuaHGItr2-HY3Oc5ofAANNWlZ6_I5LlRzY9ilil_A/pub?gid=0&single=true&output=csv'
data = pd.read_csv(url)
time  = data.iloc[1:, 1]
rate  = data.iloc[1:, 2]
total = data.iloc[1:, 4]


## Section 3 - Fit  a linear regression line
#rate_log  = np.log(rate.astype(float))            # Transform data with log
total_log = np.log(total.astype(float))
const = sm.add_constant(time)                    # Fit the model y=mx+b, not y=mx
model = sm.OLS(total_log[:-8], const[:-8]).fit()  # Negative numbers exclude the last few

# Plot the data
plt.close('all')
plt.figure(1)
plt.yscale('log')

#plt.plot(time, rate, 'o', color='tab:blue')

plt.title('Radioactivity (Low initial concentration)')
plt.xlabel('Time [s]')
plt.ylabel('Rate [counts per 13s]')


## Section 4 - Extract parameters of interest
slope     = model.params[1]
intercept = model.params[0]
uncertainty = model.conf_int(alpha=0.05)

print(slope)
print(intercept)
print(uncertainty)


## Section 5 - Create the best fit line and prediction intervals
model_prediction = model.get_prediction(const[:-8]).summary_frame(alpha=.05)
model_fit = np.exp(model_prediction['mean'])
prediction_lo = np.exp(model_prediction['obs_ci_lower'])
prediction_hi = np.exp(model_prediction['obs_ci_upper'])


## Section 6 - Plot the data
plt.plot(time[:-8], model_fit,           color='tab:green')
plt.plot(time[:-8], prediction_lo, '--', color='tab:gray')
plt.plot(time[:-8], prediction_hi, '--', color='tab:gray')
plt.plot(time,      total,       'o',     color='tab:blue' )

plt.xlim(0, 1050)
plt.ylim(bottom = 1e0)