<jupyter_start><jupyter_text>1. Poisson Probability<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import chisquare
from scipy import stats
from scipy.special import sph_harm as Y
from scipy.special import eval_genlaguerre as L
import math
import random
#fifty random values between 0 and 1
def rnge():
    range_ = []
    i = 0
    while(i<50):
        random_decimal = random.uniform(0,1)
        range_.append(random_decimal)
        i +=1
    return range_

#Divide into bins and plot histogram
hist = np.histogram(rnge(), 10)
plt.figure(1)
plt.hist(rnge(),10)
plt.title('50 binned events between 0 and 1')

#average all the random values to find
av = sum(rnge()) / 50
print('Average value = ', av)

#plot a poisson distibution and subtract it from 1
def poissprob(N):
    k = []
    i = 0
    while(i<10):
        av = sum(rnge()) / 50
        x = ((av**N)*np.exp(-av))/(math.factorial(N))
        k.append(x)
        i += 1
    return x
    
plt.figure(2)
N = 1
print(poissprob(N))
plt.plot(N, poissprob(N))<jupyter_output><empty_output><jupyter_text>2. Central Limit theorem<jupyter_code># A) 10 random values and the average of those values
def rnge(n):
    range_ = []
    i = 0
    while(i<n):
        random_decimal = random.uniform(-.5,.5)
        range_.append(random_decimal)
        i += 1
    return range_

n = 10
av = sum(rnge(n))/n
print('Average of ten values =', av)
print('-----------')



# B) 10 averages their mean and variance
def ave(n,nav):
    k = []
    i = 0
    while(i<nav):
        av = sum(rnge(n)) / n
        k.append(av)
        i +=1
    return(k)
    
n, nav = 10, 10
k = ave(n,nav)
print('Array of 10 averages = ', k)
print('-----------')
   
    #mean and variance of averages

def mean(nav):
    k = sum(ave(n,nav))/nav
    return(k)

print('MEAN OF 10 AVERAGES =', mean(nav))
print('-----------')

def var(nav):
    k = []
    i = 0
    p = ave(n,nav)
    while(i<nav):
        var = (p[i] * mean(nav))**2
        k.append(var)
        i += 1
        return(k)
    s = (1/(nav-1)) * sum(k)
    return(s)

print('VARIANCE OF 10 AVERAGES = ',var(nav))
print('-----------')
print('''Yes the variance value is small and it is what I expected.
The averages of 10 random numbers from the same set should be very close to eachother. 
So the variance between the average should be very samll.''')
print('-----------')


#C) 100 averages
n, nav = 10, 100
print('MEAN OF 100 AVERAGES =', mean(nav))
print('-----------')
print('VARIACE OF 100 AVERAGES =', var(nav))
print('-----------')
print('''The mean got smaller as the averages are more evenly distributed above and below 0.
But because there are more values the distribution widens, which is why the variance increases.
This is what I expected.''')



#D) Gaussian
n, nav = 10, 100
k = ave(n,nav)
d = mean(nav)
c = var(nav)[0]

x = np.linspace(d - 3*c, d + 3*c, 100)
a = 1/(c*np.sqrt(2*np.pi))

plt.figure(1)
plt.hist(k,35)
plt.title('35 binned events between -0.5 and 0.5')
print('This distribution roughly approximates a gaussian')

plt.figure(2)
plt.title('Gaussian Distribution')
plt.plot(x, stats.norm.pdf(x, d, c))


#plt.plot(x,Gauss(x,a,d,c), 'b-')
#plt.ylim(0,a)<jupyter_output><empty_output><jupyter_text>3. Hypothesis testing<jupyter_code>print('''A z test quantifies how far a single value is from the average value of the distribution.''')
print('----------')
print('''Chisquared quantifies how far a set of data points are from a given model, or put another way,
how well the model fits the data. A chisquared uses the z value in its calculations, but instead of 
quantifying a single value, it quantifies a function''')
print('---------')
print('11.2 is within 3 standard deviations from 8.5')
print('''if i had time i would wright a Gaussian function, plug in values and integrate from 11.2 to infinity to find the
total probablity of getting the 11.2 in a z score or higher''')<jupyter_output><empty_output><jupyter_text>4. Python Statistics<jupyter_code>#data for daily changes in price delta_P
month1_dp =[0.04, 0.04, 0.093, -0.02, 0.106, 0.013, -0.05, 0.026, 0.0, -0.01, 0.12, 0.026, 0.053, 0.093, 0.04, 0.0, 0.106, 0.186, 0.026, 0.026, -0.02, 0.053, -0.04, 0.093, 0.026, 0.053, 0.013, 0.106, -0.02, -0.04, 0.133]
month2_dp = [-0.08, 0.066, 0.0, -0.01, -0.01, 0.026, 0.013, 0.133, -0.08, -0.01, 0.12, 0.093, 0.013, -0.1, 0.026, 0.133, 0.08, 0.04, 0.013, 0.106, 0.026, 0.08, 0.08, 0.08, -0.2, 0.066, 0.106, -0.01, 0.013, 0.04, 0.026]

#making a hoistogram from data
num_bins = 12
freqs_m1, BinEdges_m1 = np.histogram(month1_dp, num_bins)
freqs_m2, BinEdges_m2 = np.histogram(month2_dp, num_bins)

bin_width_1 = BinEdges_m1[1] - BinEdges_m1[0]
x_m1 = BinEdges_m1[1:] - bin_width_1/2.0

bin_width_2 = BinEdges_m2[1] - BinEdges_m2[0]
x_m2 = BinEdges_m2[1:] - bin_width_2/2.0

#Fitting a gaussian
def mean(x,freqs):
    return(x*freqs) / sum(freqs)

def sigma(x,freqs):
    return np.sqrt(sum(freqs * x - mean(x,freqs)**2) / sum(freqs))

def Gauss(x,a,x0,sigma):
    return(a * np.exp(-(x-x0)**2 / 2 * sigma**2))

mean1 = mean(x_m1, freqs_m1)
sigma1 = sigma(x_m1,freqs_m1)
mean2 = mean(x_m2, freqs_m2)
sigma2 = sigma(x_m2, freqs_m2)

popt1, pcov1 = curve_fit(Gauss, x_m1, freqs_m1, p0 = [max(freqs_m1), mean1, sigma1], method = 'lm')
popt2, pcov2 = curve_fit(Gauss, x_m2, freqs_m2, p0 = [max(freqs_m2), mean2, sigma2], method = 'lm')

y1 = Gauss(x_m1, *popt1)
y2 = Gauss(x_m2, *popt2)

#do a reducedchisquared test
chi2_m1, pval_m1 = chisquare(freqs_m1, f_exp=y1, ddof=3)
chi2_m2, pval_m2 = chisquare(freqs_m2, f_exp=y2, ddof=3)
print(chi2_m1,pval_m1)
print(chi2_m2,pval_m2)

#Plotting
plt.style.use('seaborn-deep')
#set up xaxis bins to lump counts into eqiual sized buckets
num_bins_fig = 15
fig_bins = np.linspace(-0.25,0.25, num_bins_fig)

#set plot for figure 1
plt.figure(1)
plt.xlabel('Single day price movements [$]', fontsize = 11)
plt.xlabel('Counts', fontsize = 11)
plt.hist([month1_dp, month2_dp], bins = fig_bins, label = ['Month 1', 'Month 2'])
plt.legensd(loc= 'upper right')

#set plot for figure 2
plt.figure(2)
plt.subplot(1,2,1)
plt.gca().set_title(r'$\chi^2 /\nu =$' + str(chi2_m2))
plt.xlabel('Month 2 price Movements [$]', fontsize = 11)
plt.ylabel('Frequency / Number of observations', fontsize = 11)

plt.hist(month2_dp, bins=x_m2, color = "lightcoral")
plt.plot(x_m2, y2, color = 'green', label = 'fit')

#next subplot, includes error
plt.subplot(1,2,2)
plt.gca().set_title(r'$\chi^2 /\nu =$' + str(chi2_m1))
plt.xlabel('Month 1 price Movements [$]', fontsize = 11)
plt.ylabel('Frequency / Number of observations', fontsize = 11)

plt.hist(month2_dp, bins=x_m2, color = "lightcoral")
plt.plot(x_m2, y2, color = 'green', label = 'fit')
plt.errorbar(x_m1,y1,yerr=my_y_err,ecolor ='black')<jupyter_output><empty_output><jupyter_text>5. Visualizing Electron Orbitals<jupyter_code>def hydrogen(r,theta_pol, phi_azi, n,l,m,a=1):
    f1 = (2./(n*a))**3
    f2 = math.factorial(n-l-1) / (2.*n*math.factorial(n+l))
    factor1 = np.sqrt(f1*f2)
    
    f1_r = np.exp(-r/2.)*r**l
    f2_r = L(n-l-1, 2*l + 1, r)
    factor_r = f1_r*f2_r
    
    factor_angles = Y(m,l,phi_azi,theta_pol)
    
    return factor1*factor_r*factor_angles

#convert from xyz to rthetapi
def Spherical_np(x,y,z):
    xy = x**2 + y**2
    r = np.sqrt(xy + z**2)
    theta = np.arctan2(np.sqrt(xy), z)
    phi = np.arctan2(y, x)
    return r,theta,phi

#make a meshgrid representing a 3d volume
x0bound, x1bound = -12.0, 12.0
x_res = 100
y_res = 8 #resolution
z_res = 100
x_ = np.linspace(x0bound, x1bound, x_res)
y_ = np.linspace(x0bound, x1bound, y_res)
z_ = np.linspace(x0bound, x1bound, z_res)
x, y, z = np.meshgrid(x_, y_, z_, indexing = 'ij')

#throws up an AssertionError conditions = False
assert np.all(x[:,0,0] == x_)
assert np.all(y[0,:,0] == y_)
assert np.all(z[0,0,:] == z_)

r, theta_pol, phi_azi = Spherical_np(x,y,z)

n,l,m = 4,3,1
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)

pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)

def plot_xy_crossection(pdf,z_ind):
    plt.imshow(pdf[:,:,z_ind])
    plt.show
    
def plot_xz_crossection(pdf, y_ind, cmap0 = None):
    plt.imshow(pdf[:,y_ind,:], cmap = cmap0)
    plt.show()
    
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")<jupyter_output><empty_output><jupyter_text>A)<jupyter_code>#1 
n,l,m = 2,1,1
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 0')

#2
n,l,m = 3,2,1
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 2')

#3
n,l,m = 6,4,4
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 0')

#4
n,l,m = 5,2,1
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 2')

#5
n,l,m = 7,2,1
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 2')

#6
n,l,m = 7,3,3
wvf = hydrogen(r, theta_pol, phi_azi, n, l, m, a=10.**-10)
pdf_ = np.conjugate(wvf)*wvf
pdf = pdf_.astype(float)
plot_xz_crossection(pdf,int(len(y_)/2.0), cmap0 = "cool")
print('Azimuthal stripes = 2, Radial stripes = 2')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17398820418954647, 'fraction_numerical': 0.062131380923327234, 'mean_word_length': 3.900348779272546, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 25, 'https://': 0, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '546152', 'license': 'MIT', 'n_tokens_mistral': 4660, 'n_tokens_neox': 4194, 'n_words': 1036, 'path': 'hw 3.ipynb', 'repo_name': 'wjryan/phys_133', 'size': '9834'}