## PHYS 310, HW #2, Spring 2018¶

In [1]:
import scipy as sp
from scipy import stats

import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in the  notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib

# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended,
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')

# M.L. modifications of matplotlib defaults using syntax of v.2.0
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size


### H&H Problem 3.2¶

#### (i) Normalization¶

The p.d.f. of the uniform distribution is $$P_U(x; \bar x, a) = \left\{\begin{array}[cl]\\ 1/a & \mbox{if \quad \bar x -a/2 \leq x \leq \bar x + a/2} \\ 0 & \mbox{otherwise} \end{array}\right. .$$ To show that this is normalized, I integrate it over all values of $x$: \begin{eqnarray} \int_{-\infty}^\infty PU(x)\, dx &=& \int{\bar x - a/2}^{\bar x + a/2} \frac{1}{a}\, dx \ &=& \left\vert \frac{x}{a}\right\vert_{\bar x -\frac{a} {2}}^{\bar x + \frac{a}{2}} \ &=& \frac{1}{a}\left[\left(\bar x + \frac{a}{2}\right) - \left(\bar x + \frac{a}{2}\right)\right] \ &=& \frac{1}{a}\left[\frac{a}{2} + \frac{a}{2}\right] &=& 1 \end{eqnarray} It's normalized.

#### (ii) Mean¶

The mean is \begin{eqnarray*} \bar x &\equiv& \int_{-\infty}^\infty x PU(x)\, dx \ &=& \int{\bar x - a/2}^{\bar x + a/2} \frac{x}{a}\, dx \ &=& \frac{1}{a}\left\vert \frac{x^2}{2}\right\vert_{\bar x - \frac{a}{2}}^{\bar x + \frac{a}{2}}\ &=& \frac{1}{2a}\left[\left(\bar x + \frac{a}{2}\right)^2

        - \left(\bar x - \frac{a}{2}\right)^2 \right]\\
&=& \frac{1}{2a}\left[ \bar x^2  + \bar x a + \frac{a^2}{4} - \bar x^2 + \bar x a
-\frac{a^2}{4}\right]  \\
&=& \bar x


\end{eqnarray*}

#### (iii) standard deviation¶

The variance, or square of the standard deviation is the second moment of the p.d.f. about its mean, which is simplified in the second line using \begin{eqnarray*} \bar \sigma^2 &\equiv& \int_{-\infty}^\infty (x - \bar x)^2 PU(x)\, dx \ &=& \int{\bar x - a/2}^{\bar x + a/2} \frac{1}{a} (x -\bar x)^2\, dx \ &=& \frac{1}{a}\int{\bar x - a/2}^{\bar x + a/2} \frac{1}{a} (x^2 -2x\bar x + \bar x^2)\,dx\ &=& \frac{1}{a}\left\vert \frac{x^3}{3} -x^2\bar x + x\bar x^2 \right\vert{\bar x - \frac{a}{2}}^{\bar x + \frac{a}{2}}\ &=& \frac{1}{a}\left[ \frac{1}{3}\left(\bar x + \frac{a}{2}\right)^3

                          - \frac{1}{3}\left(\bar x +  \frac{a}{2}\right)^3
+ \bar x\left(\bar x +  \frac{a}{2}\right)^2
- \bar x\left(\bar x -  \frac{a}{2}\right)^2
+ \bar x^2\left(\bar x +  \frac{a}{2}\right)
- \bar x^2\left(\bar x -  \frac{a}{2}\right) \right] \\
&=& \frac{1}{3a} \left[ \left(\bar x^3 + 3\bar x^2 \frac{a}{2}\right) \right]


\end{eqnarray*}

### H&H Problem 3.7¶

In [2]:
counts = sp.array([1,3,4,5,6,7,8,9,10,11,12,13])
occurences = sp.array([1,2,3,6,9,11,8,8,6,2,1,1])

n_expts = sp.sum(occurences)  # check
print(n_expts)

58


#### (i) Total number of counts¶

In [3]:
total = sp.dot(occurences,counts)
print("total # of counts =", total)

total # of counts = 423


#### (ii) Mean number of counts¶

This is the best estimate of the mean of the parent Poisson distribution.

In [4]:
mean = total/n_expts
print('mean counts/second =', mean)

mean counts/second = 7.29310344828


#### (i) Expected # of trials with 5 or fewer counts¶

Use the cumulative distribution function for the Poisson distribution

In [5]:
probability = sp.stats.poisson.cdf(5, mean)
print("probability of getting 5 or fewer counts =", probability)
n_expected = n_expts*probability
print("expected number of trials with 5 or fewer counts =", n_expected)

probability of getting 5 or fewer counts = 0.264848929955
expected number of trials with 5 or fewer counts = 15.3612379374


#### (ii) Expected # of trials with 20 or more counts¶

Use the cumulative distribution function to find the probability of 19 or fewer counts. The probability of 20 or more is (1 - the probability of 19 or fewer).

In [5]:
probability = 1 - sp.stats.poisson.cdf(19, mean)
print("probability of getting 20 or more counts =", probability)
n_expected = n_expts*probability
print("expected number of trials with 5 or fewer counts =", n_expected)

probability of getting 20 or more counts = 7.67424359367e-05
expected number of trials with 5 or fewer counts = 0.00445106128433


### H&H Problem 3.8¶

#### (i) Mean count rate¶

If we use the minute of our unit of time, the determination of the mean count rate is trivial: it's $R = 270\, \mbox{min}^{-1}$. If we choose seconds, it's only slightly less trivial: $R = 270/60 = 4.5\, \mbox{s}^{-1}$.

#### (ii) Error in the mean count rate¶

We are sampling from a Poisson distribution, so the error given by the standard deviation is $\sigma = \sqrt{\mu}$.

In [7]:
sigma = sp.sqrt(270)
print('sigma (min^-1) =', sigma, '; sigma (s^-1) = ', sigma/60)

sigma (min^-1) = 16.4316767252 ; sigma (s^-1) =  0.273861278753


#### (iii) Fractional error¶

In [8]:
print('fractional error =', sigma/270)

fractional error = 0.060858061945


#### (iv) Expected count in 15 minutes¶

$N = 15\, \mbox{min}\times R$

The probability of actually getting this number of counts is given by the probability mass function of the Poisson distribution.

In [9]:
n_expected = 15*270
print('expected in 15 minutes =', n_expected)

expected in 15 minutes = 4050

In [11]:
probability = sp.stats.poisson.pmf(n_expected,n_expected)
print('probability of getting exactly', n_expected,' counts =', probability)

probability of getting exactly 4050  counts = 0.00626864416478


### H&H Problem 3.9¶

In [35]:
plt.figure(1)
n = sp.linspace(0,60,61)
x = sp.linspace(0,60,301)
p = sp.stats.poisson.pmf(n,35)
g = sp.stats.norm.pdf(x,35,sp.sqrt(35))
plt.scatter(n, p, c='red', label='Poisson')
plt.plot(x, g, label='Gaussian')
plt.xlim(0,60)
plt.axhline(0)
plt.legend();

In [ ]: