<jupyter_start><jupyter_text>GeostatsPy: Basic Univariate Statistics and Distribution Plotting for Subsurface Data Analytics in Python  <NAME>, 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) Demonstration of Spatial Boostrap Contact: <NAME>, University of Texas at Austin, Geostatistics CourseSteps:1. Build an initial sample set with $ndata$ samples.  Specify locations.  Sort the data in a data only array for    an equal weighted cumulative density function (CDF) for backtransform2. Calculate the covariance matrix (between the data) for LU simulation.  This provides a fast and easy method   to calculate spaitally correlated, unconditional, Gaussian realizations at the data locations.3. Perform LU (lower upper) decomposition of the covariance matrix.  Realizations only require L x random vector.4. Draw from this initial sample set, with replacement, $ndata$ times to build a new realization of the sample.     Repeat this $nreal$ times to make realizations of the sample.  Use LU simulation to ensure spaital correlation    between samples.  Backtransform from Gaussian to original data values (use percentiles to sample).5. Calculate the statistic of interest for each realization. This demonstration considers the mean only.     We could have considered any statistic including median, 13th percentile, skew etc. 6. (and 7.) Quantify and visualize uncertainty with histograms and summary statistics.References:Efron, 1982, The jackknife, the bootstrap, and other resampling plans, Society of Industrial and Applied Math, CBMS-NSF Monographs, 3<NAME>. Resampling from stochastic simulations. Environmental and Ecological Statistics, 1:63–84, 1994. Import PackagesWe will also need some standard packages. These should have been installed with Anaconda 3.<jupyter_code># Imports 
import pprint
import math
import scipy
import scipy.linalg                # SciPy Linear Algebra Library
from scipy.stats import norm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Set the working directoryI 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).<jupyter_code>os.chdir("c:/PGE383")             # set the working directory<jupyter_output><empty_output><jupyter_text>Assume Data Values and LocationsFor demonstration I assume 4 data values in 2D and make a CDF by sorting for the back transform from Gaussian.<jupyter_code># 1. Make Spatial Dataset and Define Number of Realizations (just need a sufficient number to summarize uncertainty)
nreal = 100                                               # User Specified
# Assign Data Values - assume 4 data locations
ndata = 4                                                 # User Specified
output = np.zeros((ndata,nreal))  # x, y, value 
data = np.zeros((ndata,3))  # x, y, value 
data[0,0] = 25.0;  data[0,1] = 50.0;  data[0,2] = 2.3     # User Specified
data[1,0] = 75.0;  data[1,1] = 80.0;  data[1,2] = 5.5     # User Specified
data[2,0] = 10.0;  data[2,1] = 25.0;  data[2,2] = 7.0     # User Specified
data[3,0] = 95.0;  data[3,1] = 15.0;  data[3,2] = 8.7     # User Specified
# Set up backtransform lookup table
value = np.zeros(ndata) 
for i in range(0, ndata):
    value[i] = data[i,2]
print(data) 
value.sort()
print(value)

# data in ascending<jupyter_output>[[25.  50.   2.3]
 [75.  80.   5.5]
 [10.  25.   7. ]
 [95.  15.   8.7]]
[2.3 5.5 7.  8.7]<jupyter_text>Calculate the Covariance MatrixThe covariance matrix is calculated from the distance matrix assuming a isotropic spherical variogram with range set as $var_range$.<jupyter_code># 2. Calculate Symmetric Covariance Array - assuming variogram with spherical structure with range specified, var_range
cov = np.zeros((ndata,ndata))
var_range = 100.0                                          # User Specified
for i in range(0, ndata):
    for j in range(0, ndata):
        distance = math.sqrt(math.pow((data[i,0]-data[j,0]),2) + math.pow((data[i,1]-data[j,1]),2))
        cova = 0.0
        if(distance < var_range):
            hr = distance / var_range
            cova = 1.0 - hr*(1.5 - 0.5*hr*hr)  # spherical structure, no nugget
            cov[i,j] = cova
cov<jupyter_output><empty_output><jupyter_text>LU DecompositionPeform LU decomposition to calculate the needed lower matrix $L$.<jupyter_code># 3. Perform LU Decomposition using scipy (used tutorial at www.quantstart.com) 
P, L, U = scipy.linalg.lu(cov)
print(L);  print(U)
# Test the LU Decomposition (available for testing only)
#Test = cov - np.matmul(L,U)
#Test  # should be zeros<jupyter_output><empty_output><jupyter_text>Calculate the RealizationsThe following steps:1. Matrix multiplication of the nxn lower matrix with a 1D matrix of random standard normal values to calculate a realization with spatial correlation at data locations $i = 0,\ldots,n-1$.2. Convert to p-values with the standard normal CDF3. Transform to discrete values in the data CDF.Realizations stored in output.<jupyter_code># 4. Calculate Unconditional Realization at the Specified Locations from the Data Distribution with Replacement
#    Realizations in Gaussian space are converted to p-values and then applied to draw from the original distribution
rand = np.zeros(ndata); pvalue = np.zeros(ndata)
for ireal in range(0, nreal):    
    for i in range(0, ndata):
        rand[i] = np.random.normal()
    realization = np.matmul(L,rand)
# Convert to p-values    
    for i in range(0, ndata):
        pvalue[i] = norm.cdf(realization[i])
# Complete backtransform by selecting percentiles from original cdf (drawing with replacement from provided data values)       
        index = int(round( ndata * pvalue[i] - 0.5))
        output[i,ireal] = value[index]
idata = 0; ireal = 0
output[idata,ireal]<jupyter_output><empty_output><jupyter_text>Calculate the Target Statistic for Each Realization<jupyter_code># 5. Calculate the summary statistic of interest (mean only for this example)
# Any summary statistic may be calculated.
summary = np.zeros(nreal) 
for ireal in range(0, nreal): 
    mean = 0
    for i in range(0, ndata):
        mean = mean + output[i,ireal]    
    mean = mean / ndata
    summary[ireal] = mean
summary
# 6. Visualize the summary statistic, calculate the variance or any other measure of spread to represent uncertainty in
#    the summary statistic.
plt.hist(summary,bins=20,normed=True)
plt.title("Spatial Bootstrap Resamples of Mean")
plt.xlabel("Mean")
plt.ylabel("Frequency")
plt.show()
# 7. Summary statistics to describe uncertainty in the mean from bootstrap
columns = ['mean']
df_mean = pd.DataFrame(summary,columns=columns)
df_mean.describe()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09309141660851361, 'fraction_numerical': 0.026099092812281927, 'mean_word_length': 3.9833101529902644, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 31, 'https://': 6, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '368042', 'license': 'MIT', 'n_tokens_mistral': 2334, 'n_tokens_neox': 2143, 'n_words': 779, 'path': 'Spatial_Bootstrap.ipynb', 'repo_name': 'UtahUrbanWater/PythonNumericalDemos', 'size': '7188'}