<jupyter_start><jupyter_text>Getting StartedThis notebook is for a short demonstration of Python for statistics.The notebook covers the following points:- Generating random numbers- Probability distributions- Curve fitting- Simple statistical anlysis with NetCDF files  Import Python modules<jupyter_code># Special '%matplotlib inline' command should go first, before the matplotlib import
%matplotlib inline
import numpy as np
import scipy.stats  as stats
import matplotlib.pyplot as plt
import time
from datetime import datetime, timedelta
from pathlib import Path<jupyter_output><empty_output><jupyter_text>Random numbers<jupyter_code># Get a random integer in the specified range
np.random.randint(0, 10)  
# Get a random element from a sequence
np.random.choice([2, 4, 6, 9])
# Get a random real number between 0 and 1
np.random.random()  
# Get a random real number in the specified range
np.random.uniform(0, 10)<jupyter_output><empty_output><jupyter_text>Probability distributions https://docs.scipy.org/doc/scipy/reference/tutorial/stats/continuous.htmlhttps://docs.scipy.org/doc/numpy-1.14.0/reference/routines.random.html<jupyter_code>np.random.seed(122)   # Set the seed to an arbitrary value Set the 

[np.random.uniform(0,10) for x in range(6)]<jupyter_output><empty_output><jupyter_text>gaussian distribution with mean = 2.0, and a standard deviation of 1.0:<jupyter_code>gaussian = stats.norm(loc=4.0, scale=1.0)
x = np.linspace(0.0, 8.0, 100)
y = gaussian.pdf(x)
plt.plot(x,y);<jupyter_output><empty_output><jupyter_text>Cumulative distribution<jupyter_code>z = gaussian.cdf(x)
plt.plot(x,z);<jupyter_output><empty_output><jupyter_text>For hypothesis testing, one often needs the p-value. For example, for the given gaussian distribution above, what would be the x-value so that P(X <= x) = 0.95?<jupyter_code>gaussian.ppf(0.95)
x = gaussian.rvs(1000);
# density --> If True, the first element of the return tuple will be the counts normalized to form a probability density,
# i.e., the area (or integral) under the histogram will sum to 1.
plt.hist(x, bins=50, density=True);<jupyter_output><empty_output><jupyter_text>Curve fitting<jupyter_code># Seed the random number generator for reproducibility
np.random.seed(0)

x_data = np.linspace(-5, 5, num=50)
y_data = 2.9 * np.sin(1.5 * x_data) + np.random.normal(size=50)

# And plot it
plt.figure(figsize=(6, 4));
plt.scatter(x_data, y_data);
from scipy import optimize

def test_func(x, a, b):
    return a * np.sin(b * x)

params, params_covariance = optimize.curve_fit(test_func, x_data, y_data, p0=[2, 2])
print(params); # param[0] and param[1] are the mean and the standard deviation
plt.figure(figsize=(6, 4))
plt.scatter(x_data, y_data, label='Data')
plt.plot(x_data, test_func(x_data, params[0], params[1]),
         label='Fit function', color='r')

plt.legend(loc='best')

plt.show()
from scipy.stats import norm

# picking 500 of from a normal distrubution
# with mean 0 and standard deviation 1
samp = norm.rvs(loc=0, scale=True, size=500)  # rvs --> Random variates

param = norm.fit(samp) # distribution fitting

# generate 100 numbers between -5 and 5
x = np.linspace(-5, 5, 100)

# fitted distribution
pdf_fitted = norm.pdf(x, loc=param[0], scale=param[1])

# original distribution
pdf = norm.pdf(x)

plt.figure(figsize=(8, 6))
plt.title('Normal distribution')
plt.plot(x, pdf_fitted, 'r-', label='fitted distribution')
plt.plot(x, pdf, 'b-', label='original distribution')
plt.hist(samp, density=True, alpha=.3, label='sample')
plt.legend(loc='best', frameon=False)
plt.show()<jupyter_output><empty_output><jupyter_text>Exploring NetCDF data<jupyter_code>import netCDF4<jupyter_output><empty_output><jupyter_text>Read datasetThis uses the pathlib.Path class, which can make working with filenames, directories and paths overall slightly more convenient<jupyter_code>directory = Path("../../_data")
directory /= "cmip5/rcp85/Amon/pr/HadGEM2-ES/r1i1p1"
path = directory / "pr_Amon_HadGEM2-ES_rcp85_r1i1p1_200512-203011.nc"
dataset = netCDF4.Dataset(path, "r")<jupyter_output><empty_output><jupyter_text>Discover the dataA dataset when listed at the end of a cell already shows a lot of information.<jupyter_code>dataset<jupyter_output><empty_output><jupyter_text>Below, we access some variables separately<jupyter_code>pr = dataset['pr']
lat = dataset['lat']
lon = dataset['lon']
time = netCDF4.num2date(dataset['time'][:], dataset['time'].units)
time[0]
# Shows a tuple of a 1-element tuples, one for each shape
time.shape, lat.shape, lon.shape
# pr_Amon_HadGEM2-ES_rcp85_r1i1p1_200512-203011.nc
# 2005 December to 2030 November
pr.shape, pr.size, "Number of months = " + str((2030-2005)*12)<jupyter_output><empty_output><jupyter_text>Getting the actual dataTo access the actual data from a dataset as a NumPy array, we need to actually access and copy the data, using the `[:]` operation.All variables above, `pr`, `lat`, `lon` and `time`, are still netCDF *variables*, not arrays. This has the convenience that we can directly access their netCDF attributes, and we avoid reading the full data into memory. The disadvantage is that we can't use them fully as a NumPy array. The code below tries to show the difference.<jupyter_code>pr.ncattrs()
pr.units<jupyter_output><empty_output><jupyter_text>Trying to use it as a NumPy array will not work:<jupyter_code>try:
    pr.max()
except AttributeError as exc:
    print(exc)<jupyter_output>NetCDF: Attribute not found<jupyter_text>But, if we, or the function, *requires* the data to be an array, it does work, since `np.max` takes an array.This also forces the data to be read into memory. (Under the hood, the data will be copied somewhere.)Usually, leaving `pr` and friends as a netCDF variable works for a lot of functionality.<jupyter_code>np.max(pr)<jupyter_output><empty_output><jupyter_text>`pr` itself is still a netCDF variable<jupyter_code>pr<jupyter_output><empty_output><jupyter_text>We now *permanently* turn `pr` into a NumPy array, by forcing a copy. Note that `[:]` also works for multi-dimensional variables.  We lose all the netCDF variables attributes and methods for it, but gain the NumPy ones.<jupyter_code>pr = pr[:]<jupyter_output><empty_output><jupyter_text>`pr` is now a NumPy array. A masked array, since there may be missing values.<jupyter_code>pr<jupyter_output><empty_output><jupyter_text>There is one problem with the above `pr = pr[:]` assignment: if we nog go back a few cells in our notebook, where `pr` was still a netCDF variable, these cells won't work (apart from some functionality that NumPy arrays and netCDF variables have in common): our `pr` Python variable has changed type (and value, essentially). We would have to reassign `pr` to `dataset.variables['pr']` again for things to work normally.  So be careful with reassigning variable names in notebooks in general: it will likely make your notebook less interactive, and more (forced) linear.Let's reset `pr` to be a variable, and use that instead of the NumPy array in our further analysis<jupyter_code>pr = dataset['pr']<jupyter_output><empty_output><jupyter_text>Getting percentile values<jupyter_code>np.percentile(pr, 99), np.max(pr), np.percentile(pr, 1), np.min(pr)<jupyter_output><empty_output><jupyter_text>Plotting the data<jupyter_code>import matplotlib.pyplot as plt
import matplotlib.colors as colors
# import cartopy.crs as ccrs
# import cartopy.io.shapereader as shpreader
# from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter<jupyter_output><empty_output><jupyter_text>Find an interesting date<jupyter_code>time[178]
plt.figure(figsize=(20, 10))
plt.imshow(pr[178,:,:]);<jupyter_output><empty_output><jupyter_text>Fit 2D gaussian<jupyter_code>print(np.ndim(pr[178,:,:]))
pr[178,0:5,0:5]
# curvefit with non linear least squares (curve_fit function)
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a * np.sin(x[0]) + b * np.cos(x[1]) + c

limits = [0, 2*np.pi, 0, 2*np.pi]  # [x1_min, x1_max, x2_min, x2_max]
side_x = np.linspace(limits[0], limits[1], 100)
side_y = np.linspace(limits[2], limits[3], 100)
X1, X2 = np.meshgrid(side_x, side_y)
size = X1.shape
x1_1d = X1.reshape((1, np.prod(size)))
x2_1d = X2.reshape((1, np.prod(size)))

xdata = np.vstack((x1_1d, x2_1d))
original = (3, 1, 0.5)
z = func(xdata, *original)
Z = z.reshape(size)
z_noise = z + .2*np.random.randn(len(z))
Z_noise = z_noise.reshape(size)

ydata = z_noise
popt, pcov = curve_fit(func, xdata, ydata)
# Python 3.6 & later f-strings: insert variables directly inside a string with {variable_name}
print (f"original: {original} \nfitted: {popt}")
z_fit = func(xdata, *popt)
Z_fit = z_fit.reshape(size)

# Plot it
plt.figure(figsize=(20, 5))
plt.subplot(1, 3, 1)
plt.title("Real Function")
plt.pcolormesh(X1, X2, Z)
plt.axis(limits)
plt.colorbar()
plt.subplot(1, 3, 2)
plt.title("Function w/ Noise")
plt.pcolormesh(X1, X2, Z_noise)
plt.axis(limits)
plt.colorbar()
plt.subplot(1, 3, 3)
plt.title("Fitted Function from Noisy One")
plt.pcolormesh(X1, X2, Z_fit)
plt.axis(limits)
plt.colorbar()

plt.show()<jupyter_output>original: (3, 1, 0.5) 
fitted: [2.99607571 0.99869196 0.49719012]<jupyter_text>Extreme cases<jupyter_code>indices = pr > np.percentile(pr, 99.999)
# Note: need pr[:] here, since indices will be 3-dimensional, 
# which does not work for a netCDF variable
extreme_values = pr[:][indices]
ind_dates = np.ma.where(pr > np.percentile(pr, 99.999))[0]

extreme_dates = time[ind_dates]
plt.figure(figsize=(20, 10))
plt.plot(extreme_dates, extreme_values);<jupyter_output><empty_output><jupyter_text>Change over the years<jupyter_code>fig = plt.figure(figsize=(16, 6), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(time, pr[:, 50, 3])
plt.show();<jupyter_output><empty_output><jupyter_text>Save and close<jupyter_code>dataset.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12361379591006207, 'fraction_numerical': 0.030827144165225355, 'mean_word_length': 6.118030412744388, 'pattern_counts': {'":': 0, '<': 100, '<?xml version=': 0, '>': 103, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '494969', 'license': 'Apache-2.0', 'n_tokens_mistral': 3743, 'n_tokens_neox': 3498, 'n_words': 989, 'path': 'examples/statistics.ipynb', 'repo_name': 'RETURN-project/eucp', 'size': '9829'}