## Confidence Intervals and Hypothesis Testing in Python for Engineers and Geoscientists 
### Michael Pyrcz, Associate Professor, University of Texas at Austin 

#### Contacts: [Twitter/@GeostatsGuy](https://twitter.com/geostatsguy) | [GitHub/GeostatsGuy](https://github.com/GeostatsGuy) | [www.michaelpyrcz.com](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446)

This is a tutorial / demonstration of **Confidence Intervals and Hypothesis Testing in Python**.  In Python, the SciPy package, specifically the Stats functions (https://docs.scipy.org/doc/scipy/reference/stats.html) provide excellent tools for efficient use of statistics.  

I have previously provided these examples worked out by-hand in Excel (https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg_R.xlsx) and also in R (https://github.com/GeostatsGuy/LectureExercises/blob/master/Lecture7_CI_Hypoth_eg.R).  In all cases, I use the same dataset available as a comma delimited file (https://git.io/fxLAt).    

This tutorial includes basic, typical confidence interval and hypothesis testing methods that would commonly be required for Engineers and Geoscientists including:

1. Student-t confidence interval for the mean
2. Student-t hypothesis test for difference in means (pooled variance)
3. Student-t hypothesis test for difference in means (difference variances), Welch's t Test
3. F-distribution hypothesis test for difference in variances 

##### Caveats

I have not included all the details, specifically the test assumptions in this document.  These are included in the accompanying course notes, Lec08_hypothesis.pdf.

#### Project Goal

0. Introduction to Python in Jupyter including setting a working directory, loading data into a Pandas DataFrame.
1. Learn the basics for working with confidence intervals and hypothesis testing in Python.  
2. Demonstrate the efficiency of using Python and SciPy package for statistical analysis.

#### Load the required libraries

The following code loads the required libraries.



```python
import os                                                   # to set current working directory 
import numpy as np                                          # arrays and matrix math
import scipy.stats as st                                    # statistical methods
import pandas as pd                                         # DataFrames
```

If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  



#### Set the working directory

I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time).  Also, in this case make sure to place the required (see below) data file in this directory.  When we are done with this tutorial we will write our new dataset back to this directory.  


```python
os.chdir("C:\PGE337")                                  # set the working directory
```

#### Loading Data

Let's load the provided dataset. 'PorositySamples2Units.csv' is available at https://github.com/GeostatsGuy/GeoDataSets. It is a comma delimited file with 20 porosity measures from 2 rock units from the subsurface, porosity (as a fraction). We load it with the pandas 'read_csv' function into a data frame we called 'df' and then preview it by printing a slice and by utilizing the 'head' DataFrame member function (with a nice and clean format, see below).



```python
#df = pd.read_csv("PorositySample2Units.csv")                # read a .csv file in as a DataFrame
df = pd.read_csv(r"https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/PorositySample2Units.csv") # load data from Prof. Pyrcz's github
print(df.iloc[0:5,:])                                       # display first 4 samples in the table as a preview
df.head()                                                   # we could also use this command for a table preview 
```

         X1    X2
    0  0.21  0.20
    1  0.17  0.26
    2  0.15  0.20
    3  0.20  0.19
    4  0.19  0.13





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>X1</th>
      <th>X2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.21</td>
      <td>0.20</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.17</td>
      <td>0.26</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.15</td>
      <td>0.20</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.20</td>
      <td>0.19</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.19</td>
      <td>0.13</td>
    </tr>
  </tbody>
</table>
</div>



It is useful to review the summary statistics of our loaded DataFrame.  That can be accomplished with the 'describe' DataFrame member function.  We transpose to switch the axes for ease of visualization.


```python
df.describe().transpose()   
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>count</th>
      <th>mean</th>
      <th>std</th>
      <th>min</th>
      <th>25%</th>
      <th>50%</th>
      <th>75%</th>
      <th>max</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>X1</th>
      <td>20.0</td>
      <td>0.1645</td>
      <td>0.027810</td>
      <td>0.11</td>
      <td>0.1500</td>
      <td>0.17</td>
      <td>0.19</td>
      <td>0.21</td>
    </tr>
    <tr>
      <th>X2</th>
      <td>20.0</td>
      <td>0.2000</td>
      <td>0.045422</td>
      <td>0.11</td>
      <td>0.1675</td>
      <td>0.20</td>
      <td>0.23</td>
      <td>0.30</td>
    </tr>
  </tbody>
</table>
</div>



Here we extract the X1 and X2 unit porosity samples from the DataFrame into separate arrays called 'X1' and 'X2' for convenience.


```python
X1 = df['X1']
X2 = df['X2']
```

#### Confidence Intervals

Let's first demonstrate the calculation of the confidence interval for the sample mean at a 95% confidence level.  This could be interpreted as the interval over which there is a 95% confidence that it contains the true population.  We use the student's t distribution as we assume we do not know the variance and the sample size is small. 

\begin{equation}
x̅ \pm t_{\frac{\alpha}{2},n-1} \times \frac {s}{\sqrt{n}} 
\end{equation}


```python
ci_95_x1 = st.t.interval(0.95, len(df)-1, loc=np.mean(X1), scale=st.sem(X1))
print('The confidence interval for the X1 interval is ' + str(ci_95_x1))

```

    The confidence interval for the X1 interval is (0.1514843093952749, 0.17751569060472505)


One can check the Excel file linked above with the confidence interval calculated by hand and confirm that this result is correct.

##### Hypothesis Testing

Now, let's try the t test, hypothesis test for difference in means. This test assumes that the variances are similar along with the data being Gaussian distributed (see the course notes for more on this).  This is our test:

\begin{equation}
H_0: \mu_{X1} = \mu_{X2}
\end{equation}

\begin{equation}
H_1: \mu_{X1} \ne \mu_{X2}
\end{equation}

For the resulting t-statistic and p-value we run this command.

##### Pooled Variance t-test Difference in Means


```python
t_pooled, p_pooled = st.ttest_ind(X1,X2) # assuminng equal variance
print('The t statistic is ' + str(t_pooled) + ' and the p-value is ' + str(p_pooled))
```

    The t statistic is -2.9808897468855644 and the p-value is 0.004992130565788754


The p-value, $p$, is the symmetric interval probaiblity our outside.  In other words the $p$ reported is 2 x cumulative probability of the t statistic applied to the sampling t distribution.  Another way to look at it, if one used the $\pm t_{t_{statistic},.d.f}$ statistic as thresholds, $p$ is the probability being outside this symmetric interval. So we will reject the null hypothesis if $p \lt \alpha$.  From the p-value alone it is clear that we would reject the null hypothesis and accept the alternative hypothesis that the means are not equal.  

In case you want to compare the t-statistic to t-critical, we can apply the inverse of the student's t distribution at $\frac{\alpha}{2}$ and $1-\frac{\alpha}{2}$ to get the upper and lower critcal values.       


```python
t_critical = st.t.ppf([0.025,0.975], df=len(X1)+len(X2)-2)
print('The t crical lower and upper values are ' + str(t_critical))
```

    The t crical lower and upper values are [-2.02439416  2.02439416]


We can observe that, as expected, the t-statistic is outside the t-critcal interval.  These results are exactly what we got when we worked out the problem by hand in Excel, but so much more efficient!  

##### Welch's t-test Difference in Means 

Now let's try the t-test, hypothesis test for difference in means allowing for unequal variances, this is also known as the Welch's t test.  All we have to do is set the parameter 'equal_var' to false, note it defaults to true (e.g. the command above). 


```python
st.ttest_ind(X1, X2, equal_var = False) # allowing for difference in variance
```




    Ttest_indResult(statistic=-2.9808897468855644, pvalue=0.005502572350112331)



Once again we can see by $p$ that we will clearly reject the null hypothesis.  

##### F-test Difference in Variances

Let's now compare the variances with the F-test for difference in variances.  

\begin{equation}
H_0: \frac{\sigma^{2}_{X_2}}{\sigma^{2}_{X_1}} = 1.0
\end{equation}

\begin{equation}
H_1: \frac{\sigma^{2}_{X_2}}{\sigma^{2}_{X_1}} > 1.0
\end{equation}

Note, by ordering the variances we eliminate the case of $\sigma^{2}_{X_2} \lt \sigma^{2}_{X_1}$.

Details about the test are available in the course notes (along with assumptions such as Gaussian distributed) and this example is also worked out by hand in the linked Excel workbook.  We can accomplish the F-test in with SciPy.Stats the function with one line of code if we calculate the ratio of the sample variances ensuring that the larger variance is in the numerator and get the degrees of freedom using the len() command, ensuring that we are consistent with the numerator degrees of freedom set as 'dfn' and the denominator degrees of freedom set as 'dfd'.  We take a p-value of $1-p$ since the test is configured to be a single, right tailed test.    


```python
p_value = 1 - st.f.cdf(np.var(X2)/np.var(X1), dfn=len(X2)-1, dfd=len(X1)-1)
p_value
```




    0.01918734806315381



Once again we would clearly reject the null hypothesis since $p \lt alpha$ and assume that the variances are not equal.

#### Comments

We are just scratching the surface for confidence intervals and hypothesis tests.  Once again there are a lot of details left out of the problem formulation and assumptions, see the course notes for more coverage.  By running the same confidence interval and hypothesis tests 1) by hand in Excel and with 2) R and 3) Python code, I hope this demonstration will enable and encourage more engineers and scientists to make these R and Python tools part of their common practice. I'm always happy to discuss,

*Michael*

Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin
On twitter I'm the @GeostatsGuy.



```python

```

b1ecf2330375dd1194a5d7a16ed6e1b521863e5f
ipynb
Jupyter Notebook
PythonDataBasics_Hypothesis.ipynb
caf3676/PythonNumericalDemos
206a3d876f79e137af88b85ba98aff171e8d8e06
MIT
2017-10-15T02:07:38.000Z
2022-03-30T15:27:14.000Z
PythonDataBasics_Hypothesis.ipynb
caf3676/PythonNumericalDemos
206a3d876f79e137af88b85ba98aff171e8d8e06
MIT
2019-08-21T10:35:09.000Z
2021-02-04T04:57:13.000Z
PythonDataBasics_Hypothesis.ipynb
caf3676/PythonNumericalDemos
206a3d876f79e137af88b85ba98aff171e8d8e06
MIT
2018-06-27T11:20:30.000Z
2022-03-25T16:04:24.000Z
Qwen/Qwen-72B
1. YES
2. YES