# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ### INTRODUCTION
# ---
# * The aim of this notebook is to show how a few lines of code in python using Pandas, NumPy and Matplotlib help perform statistical analysis on a dataset with apparently minimal information.

# ### APPLICATION OF STATISTICAL INFERENCE ON A REAL DATA SET WITH AN EMPHASIS ON:
# ---
# * Poisson Process
# * Exponential Distribution
# * Bootstrap Sampling
# * Confidence Interval

# + _cell_guid="b1076dfc-b9ad-4769-8c92-a6c4dae69d19" _uuid="8f2839f25d086af736a60e9eeb907d3b93b6e0e5"
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math

plt.style.use('_classic_test_patch')
sns.set_style('whitegrid')

# + _cell_guid="79c7e3d0-c299-4dcb-8224-4455121ee9b0" _uuid="d629ff2d2480ee46fbb7e2d37f6b5fab8052498a"
data = pd.read_csv('../input/pipeline-accidents/database.csv')
cols_of_interest = ['Accident Date/Time','Accident State','Pipeline Location','Liquid Type','Net Loss (Barrels)','All Costs']
data = data[cols_of_interest]
# -

#show dataframe
data.head()

# #### SUMMARY StATISTICS OF QUANTITATIVE FEATURES

data[['Net Loss (Barrels)', 'All Costs']].describe()

# ### INSIGHTS
# * Both our variables presence of extreme values can easily be inferred because of extremely large MAX values as compared to their means.
# * Large standard deviation indicates a large variance from the mean values of the features.

#minimize the staandard deviation by dividing it to 1000000
data['All Costs'] = data['All Costs'] / 1000000
#convert data to pandas datetime object
data['Accident Date/Time'] = pd.to_datetime(data['Accident Date/Time'])

# +
fig, ax = plt.subplots(1,2,figsize=(15,5))

sns.boxplot(data['All Costs'], ax=ax[0]);
sns.boxplot(data['Net Loss (Barrels)'], ax=ax[1]);

ax[0].set_title('Accident Costs in Million USD');
ax[1].set_title('Net Loss (Barrels)');
# -

# ### INSIGHTS
# ---
# * The absence of the box in the box plot for both variables clearly indicate that the max values for both variables are significantly larger than the rest of the sample.
# * The distribution of the data is highly unreliable. With such a distribution, fitting a model to the variable is difficult unless the extreme values can be eliminated.

# ### NUMBER OF ACCIDENTS

# +
totalTimespan = data['Accident Date/Time'].max() - data['Accident Date/Time'].min()
totaltime_hour = (totalTimespan.days*24 + totalTimespan.seconds/(3600))
totaltime_month = (totalTimespan.days + totalTimespan.seconds/(3600*24)) *12/365

lmda_h = len(data) / totaltime_hour
lmda_d = len(data) / totalTimespan.days
lmda_m = len(data) / totaltime_month

print('Number of Accidents per Hour: ', lmda_h);
print('Number of Accidents per Day: ', lmda_d);
print('Number of Accidents per Month: ', lmda_m);
# -

# ### MODELLING THE FREQUENCY OF ACCIDENTS IN A CERTAIN INTERVAL

# +
PX = {}
#calculate poisson probabilities from range 0-66
for x in range(66):
    PX[x] = math.pow(2.71828, -33) * math.pow(33, x) / math.factorial(x)
p_poisson = pd.DataFrame(PX.items(), columns=['X', 'PX'])


#commulative distribution function
def cdf(data):
    n = len(data)
    x = np.sort(data)
    y = np.arange(1, n+1) / n
    return x,y

np.random.seed(101)
samples_poisson = np.random.poisson(33, 5000)
x,y = cdf(samples_poisson)

# +
#show probability mass function
fig, ax = plt.subplots(1,2, figsize=(15,5))
plt.tight_layout(3)
ax[0].set_title('PROBABILITY MASS FUNCTION', fontsize=18);
ax[0].set_xlabel('X', fontsize=14);
ax[0].set_ylabel('PX', fontsize=14);
ax[0].plot(p_poisson.X, p_poisson.PX, lw=2, marker='.', color='salmon');

#show commulative distribution function
plt.figure(figsize=(10,5))
ax[1].set_title('COMMULATIVE DISTRIBUTION FUNCTION', fontsize=18);
ax[1].set_xlabel('NUMBER OF ACCIDENTS PER MONTH', fontsize=14)
ax[1].set_ylabel('PX', fontsize=14);
ax[1].plot(x, y, lw=2, marker='.', color='steelblue');
# -

# ### INSIGHTS
# ---
# * With the probability mass function plot, we can tell that there's only a 7% chance that 33 accidents will occur in a month even though the mean number of accicents per month is 33.
# * The cumulative distribution function tells us that there is a 50% chance of 33 or less accidents occurring per month or there is a 50% chance of 33 or more accidents occurring per month which is not intuitive by only knowing the mean.

# ### MEAN TIME BETWEEN ACCIDENTS

# +
data.sort_values(by=['Accident Date/Time'], ascending = True, inplace = True)
data['timetoAccident'] = data['Accident Date/Time'].diff() #compute time between each accident
data['timetoAccident_h'] = data.apply(lambda x: x['timetoAccident'].days * 24 + x['timetoAccident'].seconds/3600, axis = 1)
data= data[data.timetoAccident.notnull()]

mean_time = np.mean(data['timetoAccident_h']) #parameter for exponential distribution
print('Mean time between accidents: ', mean_time)
# -

# ### EXPONENTIAL DISTRIBUTION
# * To model the time elapses before an event occurs.

x,y = cdf(data['timetoAccident_h'])
#get 5000 exponential smaples
samples_exp = np.random.exponential(mean_time, 5000)
x_exp, y_exp = cdf(samples_exp)

# +
plt.figure(figsize=(10,5))

plt.plot(x, y);
plt.plot(x_exp, y_exp);
plt.title('COMMULATIVE DISTRIBUTION FUCTION');
plt.xlabel('Time Between Accidents (Hour)');
plt.ylabel('PX');

plt.legend(['Actual sample', 'Exponential Dist.']);
# -

plt.figure(figsize=(10,5))
plt.title('PROBABILITY DENSITY FUNCTION');
plt.xlabel('Time Between Accidents (Hour)');
sns.distplot(samples_exp, bins=40, color='steelblue');

# ### INSIGHTS
# ---
# * The above graph shows the cumulative distributions for the sample data for accidents is very similar to the cumulative distribution of a random sample from an exponential distribution with the same mean as mean time to accidents.
# * It would be appropriate to model time between oil pipeline accidents with an exponential distribution.
# * The probability density function of oil pipeline accidents sample is an exponential decay.

# ### BOOTSTRAP SAMPLING
# * To resample the data and perform statistical inference.

# +
plt.figure(figsize = (10,5))
plt.xlim([0,100])

for i in range(100):
    samples_bs = np.random.choice(data['timetoAccident_h'], size = len(data['timetoAccident_h']))
    x_bs,y_bs = cdf(samples_bs)
    plt.plot(x_bs,y_bs, lw=1, color='salmon')

x,y = cdf(data['timetoAccident_h'])
plt.plot(x,y, lw=2, color='green',)
plt.title('COMMULATIVE DISTRIBUTION FUNCTION')
plt.xlabel('Time Between Accidents (Hour)');
plt.ylabel('PX');
plt.legend(['BOOTSTRAP SAMPLES']);


# -

# ### INSIGHTS
# ---
# * The distribution of 100 bootstrap samples is plotted on top of the actual data and as expected, the samples vary slightly from the actual sample while still closely maintaining the exponential distribution.

# ### CONFIDENCE INTERVAL

# #### BOOTSTRAP REPLICATE FUNCTION

def bs_replicate(data, func):
    samples_bs = np.random.choice(data, size=len(data))
    return func(samples_bs)
def draw_bs_replicate(data, func, size):
    replicates = np.empty(size)
    
    for i in range(size):
        replicates[i] = bs_replicate(data, func)
    return replicates


bootstrap_replicates = draw_bs_replicate(data['timetoAccident_h'], np.mean, 5000)
plt.figure(figsize=(10,5))
plt.title('BOOTSTRAP REPLICATES')
plt.xlabel('TIME BETWEEN ACCIDENTS IN HOURS')
sns.distplot(bootstrap_replicates, bins=40, color='green');

# ### INSIGHTS
# ---
# * The histogram in the above plot shows the probability density of the bootstrap replicates or the mean time to accidents when the process was repeated 5,000 times.

# ### FOR 95% CONFIDENCE INTERVAL
# * We need to find the range where 95% of times the mean of our replicates falls. That’s simply all values between the 2.5th and the 97.5th percentile of the bootstrap replicates

confidence_int = np.percentile(bootstrap_replicates, [2.5,97.5])
print('95% confidence interval: ', confidence_int[0], '-',confidence_int[1], '(Hours)')

accidents_pyear = [round((confidence_int[0]/24)*365), round((confidence_int[1]/24)*365)]
print('Number of Accidents per year: ', accidents_pyear[0], '-', accidents_pyear[1])

# ### ANNUAL COST OF OIL PIPELINE ACCIDENTS

# +
annual_cost = [(data['All Costs'].median())*321*1000000, (data['All Costs'].median())*349*1000000]

print('ANNUAL COST OF OIL PIPELINE ACCIDENTS: ', round(annual_cost[0]), '-', round(annual_cost[1]), 'USD')
# -

# ### CONCLUSION
# ---
# * With this analysis, we have 95% confidence that it takes between 21.14 and 22.93 hours for an oil pipeline accident to occur in the US or we have 95% confidence that about 321 to 349 oil pipeline accidents occur every year in the US.
# * We have 95% confidence that the annual cost incurred due to oil pipeline accidents in USA is between 7.43 to 8.06 million US dollars.

StarcoderdataJupyterScriptsDedupFiltered