<jupyter_start><jupyter_text>Lab 3 - Asking a Statistical Question  PHYS434 - Advanced Laboratory: Computational Data Analysis Professor: <NAME>  Due date: 10/23/2021 By <NAME> This week we are going to concentrate on asking a statistical question. This process almost always consists of 3+ steps: 1. Writing down in words very precisely what question you are trying to ask. 2. Translating the precise english question into a mathematical expression. This often includes determining the pdf of the background (possibly including trials), and the to integral to do to obtain a probability.  3. Coverting the probability into equivalent sigma<jupyter_code># Importing needed libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy
from scipy import stats, signal
from astropy import units as u

# This sets the size of the plot to something useful
plt.rcParams["figure.figsize"] = (15,10)
# This sets the fontsize of the x- and y-labels
fsize = 30
lsize = 24<jupyter_output><empty_output><jupyter_text>Problem 1 In our first example we are looking at the temperature reading (meta-data) associated with an experiment. For the experiment to work reliably, the temperature should be at around 12 Kelvin, and if we look at the data it is mostly consistent with 12 Kelvin to within the 0.4 degree precision of the thermometry and the thermal control system (standard deviation). However, there are times when the thermal control system misbehaved and the temperature was not near 12 K, and in addition there are various glitches in the thermometry that give anomalously high and low readings (the reading does not match the real temperature). We definitely want to identify and throw out all the data when the thermal control system was not working (and the temperature was truly off from nominal). While it is possible to have an error in the thermometry such that the true temperature was fine, and we just had a wonky reading, in an abundance of caution we want to throw those values out too.<jupyter_code>d = np.append(stats.norm.rvs(loc = 12., scale = 0.4, size = 100000), [10., 10.3, 2.1, 0., 0., 15.6, 22.3, 12.7])
fig, ax = plt.subplots(1, 1)
ax.hist(d,100, density=True)
plt.tick_params(labelsize = 24)
plt.yscale('log')
ax.set_xlabel('Temperature (K)', fontsize = fsize)
ax.set_ylabel('Probability Mass', fontsize = fsize)
ax.set_title('Temperature Distribution', fontsize = fsize, fontweight = 'bold')
plt.show()<jupyter_output><empty_output><jupyter_text>A)  1. Let's play around with the data and come up with criteria for throwing out certain data points.<jupyter_code>x = np.linspace(10, 14, 1000)
d2 = stats.norm.pdf(x, loc = 12., scale = 0.4)
fig, ax = plt.subplots(1, 1)
ax.hist(d,100, density=True)
ax.plot(x, d2, linewidth = 3)
plt.tick_params(labelsize = 24)
plt.yscale('log')
ax.set_xlabel('Temperature (K)', fontsize = fsize)
ax.set_ylabel('Probability Mass', fontsize = fsize)
ax.set_title('Temperature Distribution', fontsize = fsize, fontweight = 'bold')
plt.show()
fig, ax = plt.subplots(1, 1)
ax.hist(d,100, density=True)
ax.plot(x, d2, linewidth = 3)
ax.vlines(10, 5e-6, 1e0, color='r', linestyle = '--')
ax.vlines(14, 5e-6, 1e0, color='r', linestyle = '--')
plt.tick_params(labelsize = 24)
plt.yscale('log')
ax.set_xlabel('Temperature (K)', fontsize = fsize)
ax.set_ylabel('Probability Mass', fontsize = fsize)
ax.set_title('Temperature Distribution', fontsize = fsize, fontweight = 'bold')
plt.show()<jupyter_output><empty_output><jupyter_text>Let's suggest boundaries of values 10 and 14 (+-2 on each side of the mean) to discriminate 'bad' data points - essentially, setting these as thresholds for the data to be in between.   2. If we take the survival function of 14 under our pdf we get the following probability and sigma:<jupyter_code>norm_dist = stats.norm(loc = 12., scale = 0.4)
prob = norm_dist.sf(14)
sigma = round(-stats.norm.ppf(prob, loc=0, scale=1), 4)
sigma<jupyter_output><empty_output><jupyter_text>This seems to be a good threshold for our data - so that if the data value lies beyond five sigma from the mean of the distribution, then we will throw the data point away. 5 sigma seems to be placed right outside the majority of our distribution and excludes the data point that are outliers.Then, our statistical question becomes:_Is the probability of getting the data point in our distribution smaller than $5\sigma$?_If this is the case, we will throw out the data point.  3. We now restate our question in mathematical terms. For a data point with value $V$.<jupyter_code>def exclude_data(dist, V, sigma):
    '''
    Returns True if data point should be thrown out,
            False if it should be kept.
    '''
    Vprob = dist.sf(V)
    Vsigma = -stats.norm.ppf(Vprob, loc=0, scale=1)
    if abs(Vsigma) > sigma:
        exclude = True
    else:
        exclude = False
    return exclude<jupyter_output><empty_output><jupyter_text>We run this in a loop and get:<jupyter_code>included_array = []
excluded_array = []
for item in d:
    if exclude_data(norm_dist, item, sigma):
        excluded_array.append(item)
    else:
        included_array.append(item)
print(f'Excluded: {excluded_array}')<jupyter_output>Excluded: [2.1, 0.0, 0.0, 15.6, 22.3]<jupyter_text>4. Reminder: Our 'bad' data points are {10., 10.3, 2.1, 0., 0., 15.6, 22.3, 12.7}<jupyter_code>bad_data = [10., 10.3, 2.1, 0., 0., 15.6, 22.3, 12.7]
kept_bad_data = []
for i in included_array:
    for k in bad_data:
        if i == k:
            kept_bad_data.append(i)
bad_data, kept_bad_data
len(d) - len(bad_data), len(bad_data)
len(included_array), len(excluded_array), len(excluded_array) + len(included_array), len(kept_bad_data)<jupyter_output><empty_output><jupyter_text>We construct a truth table showing our results from above:|  | **True T** | **Bad T** || --- | --- | --- || Test Include | 100000 | 3 || Test Exclude | 0 | 5 || Total | 100000 | 8 |  B) Now, we evaluate how the omissions (throwing out 'good' data) depends on the threshold (sigma) I chose above.Since the test does not omit any good data for my threshold of $ 5\sigma $, it does not depend on the threshold if the threshold increases (to a larger sigma). However, if we decreased the threshold so that the sigma would converge onto the actual good data points - so that the width of the statistical "inclusion" is narrower than the distribution of the normal distribution in the background - and then the test would start excluding good data points.  C) There are still some data points that are 'bad' data that gets into my final distribution even after the statistical test. These are located among the distribution of the background, and so they are included as my test does not omit them since they are within the 'inclusion zone' defined by my threshold of $ \pm \: 5\sigma $.There is no way I can change my threshold - effectively the width of the inclusion zone - so that the test would not also exclude good data.  Problem 2 In this example we will be looking for asteroids. If we look at the alignment of stars on subsequent images, they don't perfectly align due to atmospheric and instrumental effects (even ignoring proper motion). The resulting distribution is two-dimensional, and for this lab let's assume it is a 2D Gaussian with 1 arcsecond RMS. Or said another way, if I histogram how far all the (stationary) stars appear to have moved I get something like:<jupyter_code>a = np.vstack((stats.norm.rvs( scale = 1, size = 100000), stats.norm.rvs( scale = 1, size = 100000)))
a.shape
fig, ax = plt.subplots(1, 1)
h = ax.hist2d(a[0,:],a[1,:],bins=100, density=True);
ax.set_aspect('equal', 'box')
plt.xlim([-3 , 3])
plt.ylim([-3 , 3])
plt.title("2D Histogram of positional uncertainty", fontsize = 24)
plt.ylabel("$\Delta$y arcseconds", fontsize = 18)
plt.xlabel("$\Delta$x arcseconds", fontsize = 18)
plt.colorbar(h[3], ax=ax);<jupyter_output><empty_output><jupyter_text>If I have a potential asteroid, it will have some true movement between the images. We would like a '5 sigma' detection of movement. What is that distance in arcseconds?  1. We know that our 2D Gaussian is related to a Rayleigh distribution, such that our Rayleigh distribution will have a standard deviation of $ \sqrt{\sigma} $ if the Gaussians have a standard deviation of $ \sigma $. Let's state our statistical question in words:_What is the distance in arcseconds that when integrated from the right up to that value (distance) corresponds to a probability of 5 'sigma'?_  2. For a value V, Rayleigh distribution of $ R(x) $ and standard normal distribution of $ N(x) $:$$ \int_{V}^{\infty}{ R(x) dx} = \int_{5\sigma}^{\infty}{ N(x) dx} $$Then, we take the $isf()$ of the first integral to find the value of V.(Thus, essentially, our mathematical question asks what is the value of V that makes this equation true.)  3.<jupyter_code>prob_5sigma = 1/(3.5e6)
sigma_gaussian = 1
sigma_rayleigh = np.sqrt(sigma_gaussian)
rayleigh = stats.rayleigh(scale = sigma_rayleigh)
det = rayleigh.isf(prob_5sigma)
det
# stats.norm.isf(rayleigh.sf(prob_5sigma))
print(f'This means that the detection of movement of 5 sigma corresponds to {det} arcseconds')<jupyter_output>This means that the detection of movement of 5 sigma corresponds to 5.489676406940512 arcseconds<jupyter_text>Problem 3 As we discussed in class, one of the key backgrounds for gamma-ray telescopes are cosmic rays. Cosmic rays are charged particles—usually protons or electrons but can include atomic nuclei such a alpha particles (helium) or iron. Because of their charge cosmic rays spiral in the magnetic field of the galaxy. From the perspective of the Earth they appear to be coming uniformly from all directions like a high energy gas, and the direction the cosmic ray is travelling when it reaches the Earth tells us nothing about where it came from because we don't know what tortured path it has taken through the galaxy to reach us. However, at trillion electron volt energies and above, the spiral loops are fairly big and the sun and the moon will block cosmic rays. This means the sun and the moon appear as holes in the cosmic ray sky (cosmic rays from that direction are absorbed).Assume in a moon sized patch on the sky we normally have a cosmic ray rate of 1 cosmic ray per minute (arrivals are random in time). We observe where the moon is for 8 hours per night (not too close to the horizon) and we observe for 15 days and see 6800 cosmic rays. Let's find the signficance of our moon shadow detection.  1. We assume the cosmic rays to follow a Poisson distribution, since we are dealing with rates of events (from cosmic rays). In this problem, we are not dealing with trials since there is no look-elsewhere effect - we are not looking for the 'brightest' candidate of our signals. Rather, we are adding our exposures together to extend the time we are observing. Thus, we are convolving our distribution over 7200 times (see cell below). However, we know that a Poisson distribution convolved with another Poisson distribution is a Poisson distribution with a mean equal to the sum of the means of the previous distributions.<jupyter_code>(8 * u.hour * 15).to(u.min)/u.min # 8 hours and 15 days<jupyter_output><empty_output><jupyter_text>We state our statistical question:_What is the probability that the "normally" occurring cosmic ray background - a Poisson distribution with mean 7200 - produces a signal of 6800 cosmic rays?_  2. We will let $S = 6800$. We start by showing the background for an 8 hour exposure (1 night).<jupyter_code>N = 7200
trials = 1
mu = 1
resolution = 1

background = stats.poisson(mu*N)

xmin, xmax = (6000, 8000)
x = np.arange(xmin, xmax+1, resolution)
cx = np.arange(xmin, xmax+1, resolution/N)
# cxstairs = (np.arange(xmin, xmax+1+0.5*resolution/N, resolution) - 0.5*resolution/N)/N
cxstairs = (np.arange(xmin, xmax+1+0.5*resolution, resolution) - 0.5*resolution)
fig, ax = plt.subplots(1, 1)
plt.tick_params(labelsize = lsize/2)
ax.stairs(background.pmf(x), cxstairs, fill=True)
ax.set_xlim([6500, 7900])
ax.set_xlabel('N cosmic rays', fontsize = fsize)
ax.set_ylabel('Probability Mass', fontsize = fsize)
ax.set_title('15 days of 8-hour exposures', fontsize = fsize, fontweight = 'bold')
plt.show()
fig, ax = plt.subplots(1, 1)
plt.tick_params(labelsize = lsize/2)
ax.stairs(background.pmf(x), cxstairs, fill=True)
ax.set_xlim([6500, 7900])
ax.set_ylim([1e-21, 1e-2])
ax.set_xlabel('N cosmic rays', fontsize = fsize)
ax.set_ylabel('Probability Mass', fontsize = fsize)
ax.set_title('15 days of 8-hour exposures', fontsize = fsize, fontweight = 'bold')
ax.set_yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>Which is the $pdf()$ of the background.<jupyter_code>Y = 6800
Y<jupyter_output><empty_output><jupyter_text>Let's describe the integral that we need to do for a 6800 cosmic ray detection.Since this value is smaller than the mean of the distribution $\mu$, we need to integrate from the left ($-\infty$) up to our value $Y = 6800$. Our integral equation then becomes:$$ \int_{-\infty}^{Y}{ P(x) dx} = \int_{\sigma}^{\infty}{ N(x) dx} $$<jupyter_code>prob_moon = (background.cdf(Y)) # We have to integrate from the left, since we are observing a deviation from the normal **less than** the mean
prob_moon
print(f'The probability of detecting 6800 cosmic rays in our observation is {prob_moon:.2e}.')<jupyter_output>The probability of detecting 6800 cosmic rays in our observation is 1.01e-06.<jupyter_text>3.<jupyter_code>sigma_moon = abs(stats.norm.ppf(prob_moon))
print(f'The sigma of our detection is {sigma_moon:.3}.')<jupyter_output>The sigma of our detection is 4.75.
{'fraction_non_alphanumeric': 0.09358955552295731, 'fraction_numerical': 0.02911838051929001, 'mean_word_length': 4.745891276864728, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 56, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '511869', 'license': 'MIT', 'n_tokens_mistral': 4460, 'n_tokens_neox': 4028, 'n_words': 1799, 'path': 'Lab3/Lab3.ipynb', 'repo_name': 'SOLHAUG13/Phys434_Autumn2021', 'size': '13648'}