<p align="center">
    

</p>

## Interactive Hypothesis Testing Demonstration

### Boostrap and Analytical Methods for Hypothesis Testing, Difference in Means

* we calculate the hypothesis test for different in means with boostrap and compare to the analytical expression

* **Welch's t-test**: we assume the features are Gaussian distributed and the variance are unequal

#### Michael Pyrcz, Associate Professor, University of Texas at Austin 

##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](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) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)

#### Hypothesis Testing

Powerful methodology for spatial data analytics:

1. extracted sample set 1 and 2, the means look different, but are they? 
2. should we suspect that the samples are in fact from 2 different populations?

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}

To test this we will calculate the t statistic with the bootstrap and analytical approaches.

#### The Welch's t-test for Difference in Means by Analytical and Empirical Methods

We work with the following test statistic, *t statistic*, from the two sample sets.

\begin{equation}
\hat{t} = \frac{\overline{x}_1 - \overline{x}_2}{\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}}
\end{equation}

where $\overline{x}_1$ and $\overline{x}_2$ are the sample means, $s^2_1$ and $s^2_2$ are the sample variances and $n_1$ and $n_2$ are the numer of samples from the two datasets.

The critical value, $t_{critical}$ is calculated by the analytical expression by:

\begin{equation}
t_{critical} = \left|t(\frac{\alpha}{2},\nu)\right|
\end{equation}

The degrees of freedom, $\nu$, is calculated as follows:

\begin{equation}
\nu = \frac{\left(\frac{1}{n_1} + \frac{\mu}{n_2}\right)^2}{\frac{1}{n_1^2(n_1-1)} + \frac{\mu^2}{n_2^2(n_2-1)}}
\end{equation}

Alternatively, the sampling distribution of the t_{statistic} and t_{critical} may be calculated empirically with bootstrap.

The workflow proceeds as:

* shift both sample sets to have the mean of the combined data, $x_1$ → $x^*_1$, $x_2$ → $x^*_2$ 

* for each bootstrap realization, $\ell=1\ldots,L$

    * perform $n_1$ Monte Carlo simulations, draws with replacement, from sample set $x^*_1$
    
    * perform $n_2$ Monte Carlo simulations, draws with replacement, from sample set $x^*_2$
    
    * calculate the t_{statistic} realization, $\hat{t}^{\ell}$ given the resulting sample means $\overline{x}^{*,\ell}_1$ and $\overline{x}^{*,\ell}_2$ and the sample variances $s^{*,2,\ell}_1$ and $s^{*,2,\ell}_2$
    
* pool the results to assemble the $t_{statistic}$ sampling distribution

* calculate the cumulative probability of the observed t_{statistic}m, $\hat{t}$, from the boostrap distribution based on $\hat{t}^{\ell}$, $\ell = 1,\ldots,L$.

Here's some prerequisite information on the boostrap.

#### Bootstrap

Uncertainty in the sample statistics
* one source of uncertainty is the paucity of data.
* do 200 or even less wells provide a precise (and accurate estimate) of the mean? standard deviation? skew? P13?

Would it be useful to know the uncertainty in these statistics due to limited sampling?
* what is the impact of uncertainty in the mean porosity e.g. 20%+/-2%?

**Bootstrap** is a method to assess the uncertainty in a sample statistic by repeated random sampling with replacement.

Assumptions
* sufficient, representative sampling, identical, idependent samples

Limitations
1. assumes the samples are representative 
2. assumes stationarity
3. only accounts for uncertainty due to too few samples, e.g. no uncertainty due to changes away from data
4. does not account for boundary of area of interest 
5. assumes the samples are independent
6. does not account for other local information sources

The Bootstrap Approach (Efron, 1982)

Statistical resampling procedure to calculate uncertainty in a calculated statistic from the data itself.
* Does this work?  Prove it to yourself, for uncertainty in the mean solution is standard error: 

\begin{equation}
\sigma^2_\overline{x} = \frac{\sigma^2_s}{n}
\end{equation}

Extremely powerful - could calculate uncertainty in any statistic!  e.g. P13, skew etc.
* Would not be possible access general uncertainty in any statistic without bootstrap.
* Advanced forms account for spatial information and sampling strategy (game theory and Journel’s spatial bootstrap (1993).

Steps: 

1. assemble a sample set, must be representative, reasonable to assume independence between samples

2. optional: build a cumulative distribution function (CDF)
    * may account for declustering weights, tail extrapolation
    * could use analogous data to support

3. For $\ell = 1, \ldots, L$ realizations, do the following:

    * For $i = \alpha, \ldots, n$ data, do the following:

        * Draw a random sample with replacement from the sample set or Monte Carlo simulate from the CDF (if available). 

6. Calculate a realization of the sammary statistic of interest from the $n$ samples, e.g. $m^\ell$, $\sigma^2_{\ell}$. Return to 3 for another realization.

7. Compile and summarize the $L$ realizations of the statistic of interest.

This is a very powerful method.  Let's try it out and compare the result to the analytical form of the confidence interval for the sample mean. 


#### Objective 

Provide an example and demonstration for:

1. interactive plotting in Jupyter Notebooks with Python packages matplotlib and ipywidgets
2. provide an intuitive hands-on example of confidence intervals and compare to statistical boostrap   

#### Getting Started

Here's the steps to get setup in Python with the GeostatsPy package:

1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). 
2. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. 

#### Load the Required Libraries

The following code loads the required libraries.


```python
%matplotlib inline
from ipywidgets import interactive                        # widgets and interactivity
from ipywidgets import widgets                            
from ipywidgets import Layout
from ipywidgets import Label
from ipywidgets import VBox, HBox
import matplotlib.pyplot as plt                           # plotting
import numpy as np                                        # working with arrays
import pandas as pd                                       # working with DataFrames
from scipy import stats                                   # statistical calculations
import random                                             # random drawing / bootstrap realizations of the data
```

#### Make a Synthetic Dataset

This is an interactive method to:

* select a parametric distribution

* select the distribution parameters

* select the number of samples and visualize the synthetic dataset distribution


```python

# interactive calculation of the sample set (control of source parametric distribution and number of samples)
l = widgets.Text(value=' Interactive Hypothesis Testing, Difference in Means, Analytical & Bootstrap Methods, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))

n1 = widgets.IntSlider(min=0, max = 100, value = 10, step = 1, description = '$n_{1}$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
n1.style.handle_color = 'red'

m1 = widgets.FloatSlider(min=0, max = 50, value = 3, step = 1.0, description = '$\overline{x}_{1}$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
m1.style.handle_color = 'red'

s1 = widgets.FloatSlider(min=0, max = 10, value = 3, step = 0.25, description = '$s_1$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
s1.style.handle_color = 'red'

ui1 = widgets.VBox([n1,m1,s1],)                               # basic widget formatting 

n2 = widgets.IntSlider(min=0, max = 100, value = 10, step = 1, description = '$n_{2}$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
n2.style.handle_color = 'yellow'

m2 = widgets.FloatSlider(min=0, max = 50, value = 3, step = 1.0, description = '$\overline{x}_{2}$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
m2.style.handle_color = 'yellow'

s2 = widgets.FloatSlider(min=0, max = 10, value = 3, step = 0.25, description = '$s_2$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
s2.style.handle_color = 'yellow'

ui2 = widgets.VBox([n2,m2,s2],)                               # basic widget formatting 

L = widgets.IntSlider(min=10, max = 1000, value = 100, step = 1, description = '$L$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
L.style.handle_color = 'gray'

alpha = widgets.FloatSlider(min=0, max = 50, value = 3, step = 1.0, description = '$α$',orientation='horizontal',layout=Layout(width='300px', height='30px'),continuous_update=False)
alpha.style.handle_color = 'gray'

ui3 = widgets.VBox([L,alpha],)                                # basic widget formatting 

ui4 = widgets.HBox([ui1,ui2,ui3],)                               # basic widget formatting 

ui2 = widgets.VBox([l,ui4],)

def f_make(n1, m1, s1, n2, m2, s2, L, alpha):                                      # function to take parameters, make sample and plot

    
    np.random.seed(73073)
    x1 = np.random.normal(loc=m1,scale=s1,size=n1)
    np.random.seed(73074)
    x2 = np.random.normal(loc=m2,scale=s2,size=n2)
  
    mu = (s2*s2)/(s1*s1)
    nu = ((1/n1 + mu/n2)*(1/n1 + mu/n2))/(1/(n1*n1*(n1-1)) + ((mu*mu)/(n2*n2*(n2-1))))
    
    prop_values = np.linspace(-8.0,8.0,100)
    analytical_distribution = stats.t.pdf(prop_values,df = nu) 
    analytical_tcrit = stats.t.ppf(1.0-alpha*0.005,df = nu)
    
    # Analytical Method with SciPy
    t_stat_observed, p_value_analytical = stats.ttest_ind(x1,x2,equal_var=False)
    
    # Bootstrap Method
    global_average = np.average(np.concatenate([x1,x2]))   # shift the means to be equal to the globla mean
    x1s = x1 - np.average(x1) + global_average
    x2s = x2 - np.average(x2) + global_average
    
    t_stat = np.zeros(L); p_value = np.zeros(L)
    
    random.seed(73075)
    for l in range(0, L):                      # loop over realizations
        samples1 = random.choices(x1s, weights=None, cum_weights=None, k=len(x1s))
        #print(samples1)
        samples2 = random.choices(x2s, weights=None, cum_weights=None, k=len(x2s))
        #print(samples2)
        t_stat[l], p_value[l] = stats.ttest_ind(samples1,samples2,equal_var=False)
      
    bootstrap_lower = np.percentile(t_stat,alpha * 0.5)
    bootstrap_upper = np.percentile(t_stat,100.0 - alpha * 0.5)
    
    plt.subplot(121)
    #print(t_stat)
   
    plt.hist(x1,cumulative = False, density = True, alpha=0.4,color="red",edgecolor="black", bins = np.linspace(0,50,50), label = '$x_1$')
    plt.hist(x2,cumulative = False, density = True, alpha=0.4,color="yellow",edgecolor="black", bins = np.linspace(0,50,50), label = '$x_2$')
    plt.ylim([0,0.4]); plt.xlim([0.0,30.0])
    plt.title('Sample Distributions'); plt.xlabel('Value'); plt.ylabel('Density')
    plt.legend()
    
    #plt.hist(x2)
    
    plt.subplot(122)
    plt.ylim([0,0.6]); plt.xlim([-8.0,8.0])
    plt.title('Bootstrap and Analytical $t_{statistic}$ Sampling Distributions'); plt.xlabel('$t_{statistic}$'); plt.ylabel('Density')
    plt.plot([t_stat_observed,t_stat_observed],[0.0,0.6],color = 'black',label='observed $t_{statistic}$')
    plt.plot([bootstrap_lower,bootstrap_lower],[0.0,0.6],color = 'blue',linestyle='dashed',label = 'bootstrap interval')
    plt.plot([bootstrap_upper,bootstrap_upper],[0.0,0.6],color = 'blue',linestyle='dashed')
    plt.plot(prop_values,analytical_distribution, color = 'red',label='analytical $t_{statistic}$')
    plt.hist(t_stat,cumulative = False, density = True, alpha=0.2,color="blue",edgecolor="black", bins = np.linspace(-8.0,8.0,50), label = 'bootstrap $t_{statistic}$')

    plt.fill_between(prop_values, 0, analytical_distribution, where = prop_values <= -1*analytical_tcrit, facecolor='red', interpolate=True, alpha = 0.2)
    plt.fill_between(prop_values, 0, analytical_distribution, where = prop_values >= analytical_tcrit, facecolor='red', interpolate=True, alpha = 0.2)
    ax = plt.gca()
    handles,labels = ax.get_legend_handles_labels()
    handles = [handles[0], handles[2], handles[3], handles[1]]
    labels = [labels[0], labels[2], labels[3], labels[1]]

    plt.legend(handles,labels,loc=1)
    
    
    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)
    plt.show()


# connect the function to make the samples and plot to the widgets    
interactive_plot = widgets.interactive_output(f_make, {'n1': n1, 'm1': m1, 's1': s1, 'n2': n2, 'm2': m2, 's2': s2, 'L': L, 'alpha': alpha})
interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating
```

### Boostrap and Analytical Methods for Hypothesis Testing, Difference in Means

* including the analytical and bootstrap methods for testing the difference in means
* interactive plot demonstration with ipywidget, matplotlib packages

#### Michael Pyrcz, Associate Professor, University of Texas at Austin 

##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](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) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)

### The Problem

Let's simulate bootstrap, resampling with replacement from a hat with $n_{red}$ and $n_{green}$ balls

* **$n_1$**, **$n_2$** number of samples, **$\overline{x}_1$**, **$\overline{x}_2$** means and **$s_1$**, **$s_2$** standard deviation of the 2 sample sets
* **$L$**: number of bootstrap realizations
* **$\alpha$**: alpha level


```python
display(ui2, interactive_plot)                            # display the interactive plot
```


    VBox(children=(Text(value=' Interactive Hypothesis Testing, Difference in Means, Analytical & Bootstrap Method…



    Output()


#### Observations

Some observations:

* lower dispersion and higher difference in means increases the absolute magnitude of the observed $t_{statistic}$

* the bootstrap distribution closely matches the analytical distribution if $L$ is large enough

* it is possible to use bootstrap to calculate the sampling distribution instead of relying on the theoretical express distribution, in this case the Student's t distribution. 


#### Comments

This was a demonstration of interactive hypothesis testing for the significance in difference in means aboserved between 2 sample sets in Jupyter Notebook Python with the ipywidgets and matplotlib packages. 

I have many other demonstrations on data analytics and machine learning, e.g. on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. 
  
I hope this was helpful,

*Michael*

#### The Author:

### Michael Pyrcz, Associate Professor, University of Texas at Austin 
*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*

With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. 

For more about Michael check out these links:

#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](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) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)

#### Want to Work Together?

I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.

* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! 

* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!

* I can be reached at mpyrcz@austin.utexas.edu.

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

#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](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) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)



```python

```

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