import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as st
from scipy.optimize import minimize

rng = np.random.default_rng(seed=1)

Task 0¶

Play around with the jupyter notebook from the lecture and see if you can get better accuracy results.

Solution: First, make a box diagram in order to see if there exists missing values.

Task 1: Least-Squares Method Analytically¶

Given are 10 value pairs $(k, y_k)$ for $k$ from 1 to 10, where the $y_k$ are Poisson-distributed measurements around expectation values $\mu_k = a\cdot k$ .
The parameter $a$ is to be estimated using the weighted least-squares method. To do this, consider the following $\chi^2$ functions:

\chi^2 = \sum_{k=1}^{10} w_k (y_k - \mu_k)^2 \quad\text{with}\quad w_k \in \left\{1, \frac{1}{y_k}, \frac{1}{\mu_k} \right\}

(1)

1.1¶

Analytically calculate the estimators $\hat{a}$ as a function of $y_k$ which minimize the corresponding $\chi^2$ functions. Are there special cases where the estimators are not defined?

...

1.2¶

Compute the distributions of the estimators $\hat{a}$ using Monte-Carlo simulation for $a = 1$ and $a = 100$ . Which estimator has the smallest bias and variance?

def generate(a, n=10, size=5000):
    k = np.arange(1, n + 1)  # Indices of the y_k, so [1,2,...,10]
    y = rng.poisson(
        a * k, size=(size, n)
    )  # draw size times n (=len(ks)) poisson numbers with mean a*k
    return k, y


estimators = {
    # lambda: anonymous functions
    "1\\;\\quad": lambda k, y: ____,
    "1/y_k": lambda k, y: ____,
    "1/\\mu_k": lambda k, y: ____,
}

for a in (1, 100):
    k, y = generate(a)
    plt.figure(figsize=(5.5, 4), layout="constrained")
    for label, estim in estimators.items():
        with np.errstate(divide="ignore"):
            ab = [estim(k, yb) for yb in y]
        ma = np.mean(ab)
        sa = np.std(ab)
        w, xe = np.histogram(ab, bins=20)
        w = w / np.sum(w) / np.diff(xe)
        plt.stairs(
            w,
            xe,
            label=f"$w_k = {label} \\quad \\bar{{a}} = {ma:3.2f} \\quad \\sigma(a) = {sa:.2f}$",
        )
    plt.title(f"a = {a}")
    plt.legend(loc="upper right")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 21
     19 for label, estim in estimators.items():
     20     with np.errstate(divide="ignore"):
---> 21         ab = [estim(k, yb) for yb in y]
     22     ma = np.mean(ab)
     23     sa = np.std(ab)

Cell In[2], line 11, in <lambda>(k, y)
      3     y = rng.poisson(
      4         a * k, size=(size, n)
      5     )  # draw size times n (=len(ks)) poisson numbers with mean a*k
      6     return k, y
      9 estimators = {
     10     # lambda: anonymous functions
---> 11     "1\\;\\quad": lambda k, y: ____,
     12     "1/y_k": lambda k, y: ____,
     13     "1/\\mu_k": lambda k, y: ____,
     14 }
     16 for a in (1, 100):
     17     k, y = generate(a)

NameError: name '____' is not defined

<Figure size 4125x3000 with 0 Axes>