ENEL 649 Fall 2021 Project

import numpy as np import matplotlib.pyplot as plt import math as math #import seaborn as sns #from scipy.stats import gaussian_kde #import scipy.stats as stats import scipy.special as ss

We will define some helper functions here to assist us in our journey of statistical analysis.

def uniform_rand(min_value, max_value, size): '''Returns a numpy array of random variable uniformly distributed between <min_value> and <max_value> and of size <size>.''' return np.random.uniform(low=min_value, high=max_value, size=size) def rand_sum(min_value, max_value, size, number): '''Returns the numpy array which is sum of <number> random variable arrays. Parameters: number(num): A numerical value greater than one. Returns: numpy array of <size> elements, with axis 1, consisting sum of <number> rand arrays.''' if (number < 1 or size < 1): print("The number and size needs to be positive value aborting!") return None x = uniform_rand(min_value, max_value, size) if number == 1: return x else: for i in range(number - 1): y = uniform_rand(min_value, max_value, size) x = np.add(x, y) return x # we need to count our number of samples and bin out data def weight_of_bin(data, bin_width): '''Gives out the weighted array equal to size of <data> array that is used in plotting histogram such that the area under the histogram equals to one''' return (np.ones_like(data)/len(data))/bin_width def is_in_range(data, min_val, max_val): '''Returns true if <data> is in between <min_val> <max_val>''' if (data >= min_val and data < max_val): return True else: return False def bin_data(data, bins, data_range=0): #def bin(data, bins, min_val, max_val) ''' Returns the counts of number of samples in "data" distributed between max_val and min_val by braking it into "bins" number of bins and an array of bin range and width of each bin''' #print("Data:", data) if (data_range != 0): min_val = 0 max_val = data_range else: min_val = np.amin(data) max_val = np.amax(data) temp = min_val bins_range = np.zeros(bins+1) bins_range[0] = min_val interval_size = max_val -min_val #print("Interval size:", interval_size) sampling_interval = interval_size/bins #print("Sampling interval:", sampling_interval) binned_data = np.zeros(bins) for i in range(bins): #print("interval", temp, "and", temp + sampling_interval) bins_range[i+1] = temp+sampling_interval for j in range(data.size): if (is_in_range(data[j],temp, temp+sampling_interval)): binned_data[i] += 1 #print("sample,", data[j], "goes in:", i, "s bin") else: continue #print("no. of samples in bin,", i, "is", binned_data[i]) temp += sampling_interval return binned_data, bins_range, sampling_interval

def data_plot(data, bin_range, bin_width, pdf=np.zeros(1), pdf_range=0): '''Plot the data in given bin range, normalizing it to the area of one. if PDF if given it will plot it on top of histogram otherwise only histogram is plotted''' fig, ax = plt.subplots(figsize = (6,4)) weight = weight_of_bin(data, bin_width) histd = ax.hist(data, weights=weight, bins=bin_range, alpha=0.2, label='Histrogram samples') #x = np.linspace(bin_range[0],bin_range[-1],np.size(data_y)) #ax.plot(x, data_y, color='slategrey', alpha=0.7, label='KDE') if (pdf.any != 0): ax.plot(pdf_range, pdf, color='darkorange', label='PDF') ax.set_xlabel("Bins Range") ax.set_ylabel("No. of samples (Normalized)") ax.legend()

def mean_it(data): '''Gives out the mean of one dimensional array''' return np.sum(data)/data.size def var_it(data, mean): '''Gives out the variance of one dimensional array''' size = data.size temp = [(x - mean) ** 2 for x in data] return np.sum(temp)/size

#Calculate the normal distribution for the function # Code taken from: https://www.askpython.com/python/normal-distribution def normal_dist(data, mean, variance): '''Returns a normal distribution with given mean and variance''' return (np.pi*np.sqrt(variance)) * np.exp(-0.5*((x-mean)/np.sqrt(variance))**2)

def minimum_bin_sample(h_k, e_k, bin_range, min_samples): ''' Ensure that there are minimum of asked samples in each bin it combine bins together to ensure that. it combines bin range too to represent new range over which samples are distributed''' #print(h_k) #print(e_k) # as long as number of samples are below desired: while(np.min(h_k) < min_samples): #print(np.min(h_k)) #print("we are getting out of for loop") # our function doesn't terminate if the last sample # have less samples than desired so we take case of it first while (h_k[-1]< min_samples): h_k[-2] = h_k[-1] + h_k[-2] e_k[-2] = e_k[-1] + e_k[-2] h_k = np.delete(h_k, -1) e_k = np.delete(e_k, -1) bin_range = np.delete(bin_range, -2) for i in range(h_k.size -1): if (h_k[i] < min_samples): h_k[i] = h_k[i] + h_k[i+1] e_k[i] = e_k[i] + e_k[i+1] h_k = np.delete(h_k, i+1) e_k = np.delete(e_k, i+1) bin_range = np.delete(bin_range, i+1) #for j in range(i+1, h_k.size-1): # h_k[j] = h_k[j+1] # e_k[j] = e_k[j+1] #print("H-k", h_k) #print("e_k:", e_k) break else: #print(i) continue #print(h_k) #print(e_k) return h_k, e_k, bin_range

1. Generate 100,000 samples of a uniform random variable distributed between -3 and +3. Generate a histogram of these samples, normalize to have the same area as a PDF. Plot your histogram and the theoretical PDF function together on the same figure. They should match.

Since, random variable is uniformly distributed between -3 and 3 it's PDF is equal to 1/6

rand = uniform_rand(-3, 3, 100000) PDF_range = [-3, 3, 3.01] PDF_rand = [0, 1/6, 0] fig, ax = plt.subplots(figsize = (6,4)) number_of_samples, bin_range, bin_width = bin_data(rand, 50) weight = weight_of_bin(rand, bin_width) hist_data = ax.hist(rand, bins=bin_range, weights=weight, color='lightsteelblue',alpha=0.5, label="Histogram Data") x = np.linspace(-3, 3, 50) y2 = np.ones_like(number_of_samples)*(1/6) ax.plot(x, y2, color="teal", label="PDF") ax.legend()

2. Generate 100,000 samples of the sum of 2, 6 and 100 uniform random variables distributed between 0 and 10. Create histograms of each sum, normalize to have the same area as a PDF and plot. For each distribution, choose the number of histogram bins that produce plots that clearly show the shape of the distribution.

rand2 = rand_sum(0, 10, 100000, 2) rand6 = rand_sum(0, 10, 100000, 6) rand100 = rand_sum(0, 10, 100000, 100) number_of_samples2, bin_range2, bin_width2 = bin_data(rand2, 50) weight2 = weight_of_bin(rand2, bin_width2) number_of_samples6, bin_range6, bin_width6 = bin_data(rand6, 50) weight6 = weight_of_bin(rand6, bin_width6) number_of_samples100, bin_range100, bin_width100 = bin_data(rand100, 50) weight100 = weight_of_bin(rand100, bin_width100)

data_plot(rand2, bin_range2, bin_width2)

data_plot(rand6,bin_range6, bin_width6)

data_plot(rand100, bin_range100, bin_width100)

3. Apply a Chi-squared goodness-of-fit test for the three 100,000 sample data vectors you created in Problem 2. Your test should calculate a confidence value for each random vector. You should find that the sum of 2 and 6 uniform distributions fail the test but the sum of 100 distributions passes

x = np.linspace(0, 20, 100000) mean2 = mean_it(rand2) var2 = var_it(rand2, mean2) normal_dist2 = normal_dist(x, mean2, var2) #plt.plot(x, normal_dist2)

x = np.linspace(0, 60, 10000) mean6 = mean_it(rand6) var6 = var_it(rand6, mean6) normal_dist6 = normal_dist(x, mean6, var6)*6 #plt.plot(x, normal_dist6)

x = np.linspace(0, 1000, 100000) mean100 = mean_it(rand100) var100 = var_it(rand100, mean100) normal_dist100 = normal_dist(x, mean100, var100)*100

e_k2 = bin_data(normal_dist2, 50) e_k6 = bin_data(normal_dist6, 50) e_k100 = bin_data(normal_dist100, 50)

c2 = np.sum(((number_of_samples2-normal_dist2)**2)/normal_dist2) c6 = np.sum(((number_of_samples6-normal_dist6)**2)/normal_dist100) c100 = np.sum(((number_of_samples100-normal_dist100)**2)/normal_dist100)

dof2 = len(number_of_samples100) -2 dof6 = len(number_of_samples6) -2 dof100 = len(number_of_samples100) -2

confidence2 = 1-ss.gammainc(dof2/2, c2/2) confidence6 = 1-ss.gammainc(dof6/2, c6/2) confidence100 = 1-ss.gammainc(dof100/2, c100/2)

print(f'Confidence value for rand2 : {confidence2: .2g}') print(f'Confidence value for rand6 : {confidence6: .2g}') print(f'Confidence value for rand100 : {confidence100: .2g}')

Confidence value for rand2 :  0
Confidence value for rand6 :  0
Confidence value for rand100 :  0

4. Create a 10,000 waveform ensemble of a stochastic process where the waveform is a sinusoid with unit amplitude and random phase uniformly distributed over 0 to 2π.Your ensemble should contain 10 periods of the sinusoids and use 20 samples per period.

phase = uniform_rand(0, 2*np.pi, 10000) rad = np.linspace(0, 20*np.pi, 200)

ens = np.append([rad + phase[0]], [rad + phase[1]], axis=0) for i in range(2, np.size(phase)): ens = np.append(ens, [rad + phase[i]], axis=0) waveform_ens = np.sin(ens)

#for i in range(np.shape(ens)[0]): # plt.plot(ens[i], np.sin(ens[i])) # plt.show()

5. Use your 10,000 waveform ensemble from Problem 4 to numerically calculate a histogram that represents the first order PDF of this stochastic process. Normalize your histogram to have the same area as a PDF and plot your histogram on the same figure as the theoretical expression for the first order PDF for this stochastic process. They should match.

norm_hist = np.zeros((1,200)) for i in range(np.shape(ens)[0]): norm_hist = np.add(norm_hist, ens[i]) norm_hist = norm_hist[0]/ens.shape[0]

# randomaly choose an ensmble and plot its histogram we have 10000 x = np.random.randint(10000) fx = np.sin(ens[x]) hist_bin, hist_range, hist_width = bin_data(fx, 20) x = np.linspace(-1,1,200) pdf=1/(np.pi*np.sqrt(1-x*x))

/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in true_divide

#weight_of_bin(ens[0, ]) #plt.plot(x, pdf) data_plot(fx, hist_range, hist_width, pdf, x) #plt.show()

6. Use your 10,000 waveform ensemble from Problem 4 to numerically calculate the auto correlation function of the stochastic process. Plot the numerical auto correlation along with the theoretical auto correlation expression on the same figure. They should match

we will shift our ensemble by one then add all realizations of sine waves in ensemble and average them. The other one is just the average of all ensemble.

shifted_ens = np.ones_like(ens) for i in range(ens.shape[0]): shifted_ens[i] = ens[i] +1 shifted_avg_ens = np.zeros(200) for i in range(shifted_ens.shape[0] -1): shifted_avg_ens += shifted_ens[i] + shifted_ens[i+1] shifted_avg_ens = shifted_avg_ens/10000 avg_ens = np.zeros(200) for i in range(shifted_ens.shape[0] -1): avg_ens += ens[i] + ens[i+1] avg_ens = avg_ens/10000

y1 = np.sin(avg_ens) y2 = np.sin(shifted_avg_ens) plt.title("Average of all ensemble") plt.xlabel("Time") plt.ylabel("Sine magnitude") plt.plot(avg_ens,y1) plt.show() plt.title("Average of all ensemble shifted by a constant") plt.plot(shifted_avg_ens, y2) plt.xlabel("Time") plt.ylabel("Sine magnitude") plt.show()

Since the correlation of function is just dependent of Tau, and Tau can vary in the range of T2 - T1, which for our ensemble can be in between -22*pi to 22*pi. Times cannot be negative, but T2 can occur before T1

Furthermore it can be shown from above data, that a cut in data at constant duration, but for different realizations yields same value. We only need to for time difefference, rather than all combinations of T1 and T2 (tau)

for i in range(1, 10): print(ens[:,i+4] - ens[:,i])

[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]
[1.26295182 1.26295182 1.26295182 ... 1.26295182 1.26295182 1.26295182]

So, the correlation of the function from its theoretical  distribution can be plotted by function of tau

tau_range= np.linspace(-22*np.pi, 22*np.pi, 1000) tau_func = np.sin(tau_range) plt.title("Theoratical correlation") plt.plot(tau_range, tau_func) plt.xlabel("Tau") plt.ylabel("Correlation") plt.show()

From our ensemble, lets cut through all space by first keeping T1 fixed, then we will do it for T2, it would give us all values of Tau, since as shown above difference is same across all elements of returned array, we will just take its, first element and save in in you timecut variable

timecut_T1_ls_T2 = np.zeros(200) T1 = 0 for i in range(200): timecut_T1_ls_T2[i] = np.sin(ens[:,i] - ens[:,0])[0] timecut_T2_ls_T1 = np.zeros(200) T2 = 0 for i in range(200): timecut_T2_ls_T1[i] = np.sin(ens[:,0] - ens[:,i])[0]

now we will plot our cuts in time for T2 < T1 range will be 0 to -22pi, for T1 < T2 will be 0 to 22Pi, and see if it matches with above.

T1_lt_T2 = np.linspace(0, 22*np.pi, 200) T2_lt_T1 = np.linspace(0, -22*np.pi, 200) plt.plot(T1_lt_T2, timecut_T1_ls_T2) plt.plot(T2_lt_T1, timecut_T2_ls_T1) plt.show()