<jupyter_start><jupyter_text>06 Hypothesis and inferencePart of ["Introduction to Data Science" course](https://github.com/kupav/data-sc-intro) by <NAME>, [<EMAIL>](mailto:<EMAIL>)Recommended reading for this section:1. Gr<NAME>. (2019). Data Science From Scratch: First Principles with Python (Vol. Second edition). Sebastopol, CA: O’Reilly MediaThe following Python modules will be required. Make sure that you have them installed.- `matplotlib`- `numpy`- `scipy`  Lesson 1  Discovering properties of a deterministic systemWhen we explore a deterministic system. usually we can easily establish that two its characteristics depends on each other: We need to measure and record one and another and thus recover their dependence.Consider a system that includes a heated pot of water, a bunch of fresh eggs and a timer. We want to know how to boil the eggs. ![boil.svg](attachment:boil.svg)Conducting the experiment we boil the water, put their an egg, boil for some time, pull the egg out and check it.In this experiment we have two characteristics: time of boiling and the resulting state of eggs. Since the system is deterministic we can repeat the experiment many times varying time of boiling and discover that after 3 minutes of boiling we obtain a soft-boiled eggs and 10 minutes is enough for hard-boiled eggs.  Discovering for stochastic systemsBut for a stochastic system  we can not have solid answers and predictions.Instead we can form hypothesis and test their likelihood.Assume that we want to predict a weather using a barometer.![barometer.svg](attachment:barometer.svg)Our system includes the barometer and surrounding conditions. The experiment includes recording the barometer readings and subsequent weather observation during several next ours. Repeating the experiment many times we will discover that often a low atmospheric pressure results in a bad weather and a good weather follows a high atmospheric pressure. But we also will notice that this is not the strict dependence: in some cases high pressure can be followed by a bad weather and vice versa.In the other words the weather probably depends on the atmospheric pressure but this is not one hundred percent correct.We have an uncertainty and the only thing that can be done is to use statistical methods to estimate its degree.  Hypothesis testingDealing with stochastic experiments and trying to extract as much information from them as possible we form statistical hypothesis about the presence of some properties or dependencies and test them. This part of statistics is called hypothesis testing.It includes the following steps.- Form a null hypothesis. This is the statement that the property we are interested for is absent, one characteristic is not influenced by another and so on. This hypothesis is usually denoted as $H_0$.  - Form an alternative hypothesis $H_1$ that states the presence of the property. $H_0$ and $H_1$ must be mutually exclusive (no result could satisfy both conditions) and exhaustive (all possible results will satisfy one of the two conditions). - The purpose of the hypothesis testing is to conclude weather the null hypothesis can be rejected or not. Rejecting $H_0$ often means that we have discovered some new feature and this probably not due to chance.  - Make a decision about the significance $\alpha$. Errors are unavoidable and $\alpha$ is the accepted probability to make a type 1 error. - This error means the rejecting $H_0$ though it is true ("Yay, we've got it" by mistake). - Usual value for $\alpha$ is $0.05$, however there are   studies with $\alpha=0.01$ or even smaller.- Compute the range of the experiment outcomes such that when  the observations fall within we have to accept $H_0$ hypothesis.- For the computed range estimate the type 2 error probability   $\beta$. - This error means that we fail to reject $H_0$ even when it is false. (The discovery was right in front of us but we've missed it.)- Usually taking very small $\alpha$ results in increasing $\beta$. The overall setup is correct if $\beta$ is not much higher then $\alpha$. For $\alpha=0.05$ the appropriate is $\beta\approx 0.1$.  Formulation of an exampleBelow we will consider an example: there is a coin and we want to check if it is unfair. It means that the probability of showing up its heads or tails is not $0.5$.The only think we can do is to toss it many times and count the appearances of the heads.To draw a conclusion we need to know what ranges of values confirm the fairness or unfairness of the coin. Before proceeding with the hypothesis testing we need to model our experiment and discuss the required computer tools.  Reminder: Bernoulli variable and binomial distributionLet us first remember from the previous lecture that  mathematically our coin can be represented as a Bernoulli random variable: a discrete random variable with two outcomes: "tails" (0) or "heads" (1). The probability of heads is $P$ and tails are showed up with the probability $1-P$.To reveal the unfairness we will repeat the experiment $N$ times and count heads outcomes. Let it be $x$. For infinitely large $N\to\infty$ the count of heads tends to $x\to PN$. For example for a fair coin this is $x\to N/2$. But if $N$ is not infinite this $x$ is itself a random variable. It means that if we will repeat a series of $N$ tossing again and again we will obtain different $x$.Let us also remember: The discrete random variable $x$ is governed by the binomial distribution.   Reminder: Central limit theorem and normal distributionsDue to the central limit theorem if $N$ is sufficiently large the binomial distribution can be approximately described by a normal distribution. The PDF (probability density function) for the normal distribution is$$\rho(x)=\frac{1}{\sigma\sqrt{2\pi}}\mathrm{exp}\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$and the corresponding CDF (cumulative distribution function) reads$$\phi(x)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sigma\sqrt{2}}\right)\right]$$where the mean $\mu$ and the standard deviation $\sigma$ are $$\mu=NP,\; \sigma=\sqrt{N P (1-P)}$$These formulas can be implemented as Python functions as follows:<jupyter_code>import numpy as np
from scipy.special import erf

def get_mu_sig(N, P):
    """
    mu and sig corresponding to a binomial distribution"""
    mu = P * N
    sig = np.sqrt(N * P * (1 - P))
    return mu, sig

def norm_pdf(x, mu, sig):
    """Normal probability density function"""
    return np.exp(-(x-mu)**2 / (2*sig*sig)) / (sig * np.sqrt(2*np.pi))

def norm_cdf(x, mu, sig):
    """Normal cumulative distribution function"""
    return 0.5 * (1 + erf((x-mu)/(sig*np.sqrt(2))))<jupyter_output><empty_output><jupyter_text>Below are the plots of PDF and CDF approximating the binomial distribution at $N=1000$ and $P=0.5$.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

N = 1000
P = 0.5

mu, sig = get_mu_sig(N, P)

x = np.linspace(mu - 100, mu + 100, 100)

y1_pdf = [norm_pdf(xi, mu, sig) for xi in x]
y1_cdf = [norm_cdf(xi, mu, sig) for xi in x]

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(14, 4.5))

axs[0].plot(x, y1_pdf, label=f"mu={mu:.2f}, sig={sig:.2f}")
axs[0].set_title("PDF")
axs[0].set_xlabel(r"$x$")
axs[0].set_ylabel(r"$\rho(x)$")

axs[1].plot(x, y1_cdf, label=f"mu={mu:.2f}, sig={sig:.2f}")
axs[1].set_title("CDF")
axs[1].set_xlabel(r"$X$")
axs[1].set_ylabel(r"$P(x<X)$")

for ax in axs.reshape(-1):
    ax.set_xlabel(r'$x$')
    ax.legend()
    ax.grid()<jupyter_output><empty_output><jupyter_text>Reminder: Computing probabilities of a random normal variable with CDFLet us remember: CDF gives the probability that a random variable $x$ is less than or equal to a certain value $X$. $$P(x\leq X) = \phi(X)$$![%D1%81df_example1.svg](attachment:%D1%81df_example1.svg) The complement to CDF is the probability that a random variable $x$ is above $X$:$$P(x>X) = 1 - \phi(X)$$![%D1%81df_example2.svg](attachment:%D1%81df_example2.svg) The probability to find $x$ between $X$ and $X+h$ equals$$P(x\in [X, X+h])=\phi(X+h) - \phi(X)$$![%D1%81df_example3.svg](attachment:%D1%81df_example3.svg)  Technical note: Function inverseKnowing the CDF we can compute probabilities for a given value of a random variable. But to perform the hypothesis test we take a probability of making an error and compute the corresponding values of a random variable. It means that we need to find the inverse of the CDF.Consider a simple example: Given the exponential function$$y = e^x$$find its inverse.In this particular case it is easy since we know the inverse for the exponential function:$$x = \log y$$Below are the graphs of the original and the inverted functions.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

xx = np.linspace(-2, 5, 100)
yy = [np.exp(xi) for xi in xx]
xinv = [np.log(yi) for yi in yy] 

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))

axs[0].plot(xx, yy)
axs[0].set_title(r"$y=e^x$")
axs[0].set_xlabel(r"$x$")
axs[0].set_ylabel(r"$y$")
axs[0].grid()

axs[1].plot(yy, xinv)
axs[1].set_title(r"$x=\log y$")
axs[1].set_xlabel(r"$y$")
axs[1].set_ylabel(r"$x$")
axs[1].grid()

# Put markers to improve clarity
mx = [-1, 4.5]
my = [np.exp(-1), np.exp(4.5)]
axs[0].scatter(mx, my, c=['r', 'g'], s=[100,100])
axs[1].scatter(my, mx, c=['r', 'g'], s=[100,100]);<jupyter_output><empty_output><jupyter_text>General case: Given a function$$y = f(x)$$find the inverse$$x = f^{-1}(y)$$In general case the inverse is unknown and it can be found only numerically using computer algorithms for root finding.  Technical note: Root findingRoot finding: Given $y$ and a function $f$ find such $x$ that$$f(x) = y$$The value of $x$ that fulfills the equation is called root of this equation.Root finding routines usually seek for zeros of some function $F(x)$. It varies $x$ until finds such one that $$F(x) = 0$$To find the function inverse we take $y$, feed the root finding routine with $$F(x) = f(x) - y$$and compute the corresponding $x$.There are a lot of root finding methods. We will use `root_scalar` form the module `scipy.optimize`.To see how it works let us find such $x$ that $e^x=50$. We will use the simplest methods named "bisection". It requires preliminary root bracketing, i.e. finding a range where the root is located.As one can see in the graph below the root is near $x=4$. The values 1 and 5 are good as the root brackets.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

xx = np.linspace(1, 5, 100)
yy = [np.exp(xi) for xi in xx]

y_find = 50

fig, ax = plt.subplots()

ax.plot(xx, yy)
ax.plot([1, 5], [y_find, y_find], 'k--')
ax.set_title(r"$y=e^x$")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax.grid()<jupyter_output><empty_output><jupyter_text>Thus we can compute the root as follows.<jupyter_code>from scipy.optimize import root_scalar
import numpy as np

y = 50

sol = root_scalar(lambda x: np.exp(x) - y, method='bisect', bracket=(2, 5))
print(sol)<jupyter_output><empty_output><jupyter_text>Routine `root_scalar` returns an object that contains a lot of service information. The root is stored in its attribute `.root`. Let us test that the exponential of `sol.root` is indeed `y_find`:<jupyter_code># Test the root, must be equal to y_find
print(y_find, np.exp(sol.root))<jupyter_output><empty_output><jupyter_text>Technical note: Lambda-functionsObserve that above we pass an equation $(e^x-y)$ to the root finder as an anonymous $\lambda$-function.$\lambda$-function is a fast way to build a function right in place. Here is a simple example:<jupyter_code>def plus_one(fun, x):
    """Takes a function fun and a variable x and returns fun(x) + 1
    """
    return fun(x) + 1

# Here lambda-function computes square of its argument
y = plus_one(lambda x: x**2, 10)
print(y)<jupyter_output><empty_output><jupyter_text>To build a lambda function we write a keyword `lambda` followed by a comma separated list of its parameters. Then goes a colon  and a function body after it.Lambda-function can be assigned to a variable and then it is used as an ordinary function:<jupyter_code># This function computes a sum of squares of its parameters
newfun = lambda x, y: x**2 + y**2

print(newfun(2, 3))<jupyter_output><empty_output><jupyter_text>Inverse CDF for normal distributionNow we are ready to compute the inverse of CDF of a normal distribution.<jupyter_code>from scipy.optimize import root_scalar

def inverse_norm_cdf(p, mu, sig):
    """
    Given the probability P find the corresponding upper threshold X below which 
    the normally distributed variable x falls with the probability P.
    
    If mu != 0 or sig != 1 we compute a solution for the standard CDF and then rescale it. 
    Otherwise we would need to rescale the bracketing values (-10, 10).
    """
    if mu != 0 or sig != 1:
        return mu + sig * inverse_norm_cdf(p, mu=0, sig=1)
    
    fun = lambda x: norm_cdf(x, mu=0, sig=1) - p
    sol = root_scalar(fun, method='bisect', bracket=(-10, 10))
    return sol.root

# Simple test. We know that CDF(mu) = 0.5. Thus the inversed CDF of 0.5 must return mu
mu_test, sig_test = 100, 5
print(inverse_norm_cdf(0.5, mu=mu_test, sig=sig_test))<jupyter_output><empty_output><jupyter_text>Here is the graph of CDF of a normal distribution and its inverse.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

mu, sig = 0, 1

xx = np.linspace(-10, 10, 100)
cdf = [norm_cdf(xi, mu=mu, sig=sig) for xi in xx]
inv_cdf = [inverse_norm_cdf(p, mu=mu, sig=sig) for p in cdf]

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))

axs[0].plot(xx, cdf)
axs[0].set_title("CDF")
axs[0].set_xlabel(r"$X$")
axs[0].set_ylabel(r"$P(x<X)$")

axs[1].plot(cdf, inv_cdf)
axs[1].set_title("inv CDF")
axs[1].set_xlabel(r"$P(x<X)$")
axs[1].set_ylabel(r"$X$")

for ax in axs.reshape(-1):
    ax.grid()<jupyter_output><empty_output><jupyter_text>Probability to fall in range via CDFFinally, we will need to compute a probability of a normal random variable to fall between the given bounds. Here is the illustration to remind how it is computed.![%D1%81df_example3.svg](attachment:%D1%81df_example3.svg)<jupyter_code>def norm_btw(lo, hi, mu, sig):
    """Probability for a normal random variable to fall between lo and hi.
    """
    assert lo < hi
    
    # Probability for a normal random variable to fall below lo
    p1 = norm_cdf(lo, mu, sig)
    
    # Probability for a normal random variable to fall below hi
    p2 = norm_cdf(hi, mu, sig)
    
    return p2 - p1<jupyter_output><empty_output><jupyter_text>Back to the coin tossing example: null hypothesis, alternative and significanceNow we have done all preparation and are ready to perform the hypothesis testing.Let us recall that we want to check if a coin is unfair, $P\neq 0.5$, by tossing it many times and counting a number of heads.We perform $N=1000$ tossing and $x$ is the number of heads appeared.In this case a null hypothesis $H_0$ is that $P=0.5$ (fair coin) and the alternative $H_1$ is that $P\neq 0.5$ (unfair coin).We will accept the standard value of significance $\alpha=0.05$.  The most probable range for a fair coinNow we need to compute a range for $x$ where it falls with the probability $P_\alpha=1-\alpha=0.95$ if the coin is fair.<jupyter_code>import numpy as np

p_fair = 0.5
N = 1000

# Admitted probability of a type 1 error
alpha = 0.05

# Probability for a range of values that indicates that H0 is true
p_a = 1 - alpha

# Tails - probabilty for ranges outside the area, corresponding p_a
p_tail = (1 -  p_a) / 2

mu_0, sig_0 = get_mu_sig(N=N, P=p_fair)

# Random normal variable with (mu, sig) falls below X_lo with the probability p_tail
X_lo = inverse_norm_cdf(p_tail, mu=mu_0, sig=sig_0)

# Random normal variable with (mu, sig) is above X_hi with probability p_tail
X_hi = inverse_norm_cdf(1-p_tail, mu=mu_0, sig=sig_0)
print(f"Test bounds to accept H0 ({round(X_lo)}, {round(X_hi)})")<jupyter_output><empty_output><jupyter_text>Notice that when for a given probability $p$ we compute a value as above we in fact get $p$-quantile, i.e., a point under which a certain percent of data lies. The graph below shows the PDF for a random number of heads $x$ appearing after $N=1000$ coin tossing. The filling in the middle indicates a range where $x$ falls with the probability $P_\alpha=0.95$.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

N = 1000
P = 0.5

mu, sig = get_mu_sig(N, P)

x = np.linspace(mu - 100, mu + 100, 100)
pdf = [norm_pdf(xi, mu, sig) for xi in x]

x_a = np.linspace(X_lo, X_hi, 50)
p_a = [norm_pdf(xi, mu, sig) for xi in x_a]

fig, ax = plt.subplots()

ax.fill_between(x_a, p_a, 0, color="C1", label=f"X_lo={round(X_lo)}, X_hi={round(X_hi)}")
ax.plot(x, pdf, "C0", label=f"mu={mu:.2f}, sig={sig:.2f}")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\rho(x)$")
ax.grid()
ax.legend();<jupyter_output><empty_output><jupyter_text>The routine above allows to compute the range for any given significance $\alpha$ and the corresponding $P_\alpha=1-\alpha$. For particular probability $P_\alpha=0.95$ the range can be found in a simpler way. It is known that the 95% range is $$[\mu-1.96\sigma, \mu+1.96\sigma]$$It means that$$ P(x\in[\mu-1.96\sigma, \mu+1.96\sigma]) = 0.95$$Let us check:<jupyter_code>X_lo1 = mu - 1.96 * sig
X_hi1 = mu + 1.96 * sig

print(f"Computed  95% bounds ({round(X_lo)}, {round(X_hi)})")
print(f"Estimated 95% bounds ({round(X_lo1)}, {round(X_hi1)})")<jupyter_output><empty_output><jupyter_text>To summarize, given these bounds we can conduct our experiment with tossing a coin, and if the number of heads $x$ is outside the range we reject the null hypothesis concluding that the coin is unfair. The probability that actually this is an error (type 1 error, a fake discovery) is $\alpha=0.05$.  Type 2 error probabilityIf $x$ is between `X_lo` and `X_hi` we do not reject the null hypothesis and accept that the coin is fair. What is the probability $\beta$ that this is an error (type 2 error, a missed discovery)?To be honest, the developed approach is not very good for revealing this type of errors. This is because technically the coin becomes unfair even if $P$ deviates from $0.5$ very slightly, say at $P=0.500000001$.But we probably will not notice this small deviation performing a statistical experiment.Let us assume that the coin is considered unfair only if $P\geq 0.55$ or $P\leq 0.45$.The probability $\beta$ of the type 2 error is the probability that the unfair coin will produce $x$ (number of heads after $N$ tossing) within the found range `X_lo`, `X_hi`.Let us compute this probability for $P=0.55$ and $0.45$.<jupyter_code>p_unfair1 = 0.55
mu_1, sig_1 = get_mu_sig(N=N, P=p_unfair1)

p_unfair2 = 0.45
mu_2, sig_2 = get_mu_sig(N=N, P=p_unfair2)

# Probability that a normal random variable fall within the range X_lo, X_hi
beta1 = norm_btw(X_lo, X_hi, mu_1, sig_1)
beta2 = norm_btw(X_lo, X_hi, mu_2, sig_2)

print(f"Type 2 error probaility beta({p_unfair1})={beta1:.2f}, beta({p_unfair2})={beta2:.2f}")<jupyter_output><empty_output><jupyter_text>We observe that due to the symmetry both probabilities are identical. If they were different we would take the largest one.Let us remember that above we discussed that for $\alpha=0.05$ the obtained $\beta=0.1$ is an acceptable value.Below there is a graph of PDF for $x$ (number of heads) in a series of $N=1000$ tossing of an unfair coin with $P=0.55$. The filled region shows the range \[`X_lo`, `X_hi`\]. Its area  is $\beta=0.11$.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

N = 1000

x = np.linspace(mu_1 - 100, mu_1 + 100, 100)
pdf = [norm_pdf(xi, mu_1, sig_1) for xi in x]

x_a = np.linspace(X_lo, X_hi, 50)
p_a = [norm_pdf(xi, mu_1, sig_1) for xi in x_a]

fig, ax = plt.subplots()

ax.fill_between(x_a, p_a, 0, color="C1", label=f"X_lo={round(X_lo)}, X_hi={round(X_hi)}")
ax.plot(x, pdf, "C0", label=f"mu={mu_1:.2f}, sig={sig_1:.2f}")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\rho(x)$")
ax.grid()
ax.legend();<jupyter_output><empty_output><jupyter_text>To conclude: Type 2 error will occur if the unfair coin will produce $x$ inside the range \[`X_lo`, `X_hi`\]. In this case the unfair coin will be considered as a fair one.This probability is reasonably small, $\beta=0.11$.  The powerUsually one more characteristic is used, a power of a test $1-\beta$. The meaning of power is the probability to reject a null hypothesis, i.e., the probability to make a discovery.In our case we will reveal the unfair coin with the probability $$\mathrm{Pw} = 0.89$$and it will be an error with the probability$$\alpha=0.05$$  Tuning significanceNotice that beta and power depends on the significance. When increasing by lowering $\alpha$ (because we do not want to make a fake discovery) we increase $\beta$ (i.e., magnify a chance to miss the discovery).The analysis is correct if $\alpha$ and $\beta$ do not differ much.Here is the computation for $\alpha=0.01$<jupyter_code>import numpy as np

alpha = 0.01
N = 1000

p_fair = 0.5
p_unfair = 0.55

mu_0, sig_0 = get_mu_sig(N=N, P=p_fair)
mu_1, sig_1 = get_mu_sig(N=N, P=p_unfair)

p_a = 1 - alpha
p_tail = (1 -  p_a) / 2

X_lo = inverse_norm_cdf(p_tail, mu=mu_0, sig=sig_0)
X_hi = inverse_norm_cdf(1-p_tail, mu=mu_0, sig=sig_0)
beta = norm_btw(X_lo, X_hi, mu_1, sig_1)
Pw = 1 - beta

print(f"N={N}, X_lohi=({round(X_lo)}, {round(X_hi)}), alpha={alpha}, beta={beta}, Pw={Pw}")<jupyter_output><empty_output><jupyter_text>We see that chasing for a higher significance results in unreasonably small power. Analysis with such parameters will be incorrect.Higher significance requires longer experimental series: The appropriate $\beta$ and power for $\alpha=0.01$ is obtained at $N=2000$<jupyter_code>import numpy as np

alpha = 0.01
N = 2000

p_fair = 0.5
p_unfair = 0.55

mu_0, sig_0 = get_mu_sig(N=N, P=p_fair)
mu_1, sig_1 = get_mu_sig(N=N, P=p_unfair)

p_a = 1 - alpha
p_tail = (1 -  p_a) / 2

X_lo = inverse_norm_cdf(p_tail, mu=mu_0, sig=sig_0)
X_hi = inverse_norm_cdf(1-p_tail, mu=mu_0, sig=sig_0)
beta = norm_btw(X_lo, X_hi, mu_1, sig_1)
Pw = 1 - beta

print(f"N={N}, X_lohi=({round(X_lo)}, {round(X_hi)}), alpha={alpha}, beta={beta}, Pw={Pw}")<jupyter_output><empty_output><jupyter_text>The experiment itself with fair and unfair coinsLet us now model the discussed experiment explicitly and see how often errors will occur.<jupyter_code>import numpy as np

class Coin:
    """ Model of a coin tossing. P is the probability of heads.
    """
    def __init__(self, P):
        self.P = P
        self.rng = np.random.default_rng()
        
    def toss(self):
        return 1 if self.rng.random() < self.P else 0
# An experiment with a fair coin

N = 1000
X_lo, X_hi = 469, 531 # The H0 range for alpha=0.05, N=1000

# The coin is fair
coin = Coin(0.5)
seq = [coin.toss() for _ in range(N)]
x = sum(seq)

if X_lo <= x <= X_hi:
    print(f"x={x} within the range ({X_lo}, {X_hi}), accept that the coin is fair, no error")
else:
    print(f"x={x} outside the range ({X_lo}, {X_hi}), accept that coin is unfair, type 1 error")
# An experiment with an unfair coin

N = 1000
X_lo, X_hi = 469, 531 # The H0 range for alpha=0.05, N=1000

# The coin is unfair
coin = Coin(0.55)
seq = [coin.toss() for _ in range(N)]
x = sum(seq)

if X_lo <= x <= X_hi:
    print(f"x={x} within the range ({X_lo}, {X_hi}), accept that the coin is fair, type 2 error")
else:
    print(f"x={x} outside the range ({X_lo}, {X_hi}), accept that coin is unfair, no error")<jupyter_output><empty_output><jupyter_text>Exercises1\. Using Python functions defined above compute probability that a random variable obeyed the standard normal distribution fits the range \[-1, 1\].2\. Using Python functions defined above compute such $X$ that a random variable $x$ obeyed the standard normal distribution is $x<X$ with the probability $0.25$.3\. In a series $N=1500$ of tossing of an unknown coin the heads appeared 783 times. Can this coin be treated as unfair? Answer this question performing the hypothesis testing for significance $\alpha=0.05$. Compute probabilities of type 1 and 2 errors and the power.4\. For a series of $N=1000$ tossing of a coin develop a one sided test: accept as a null hypothesis that a coin is not biased toward heads, $P\leq 0.5$, and use significance $\alpha=0.05$  Lesson 2  $p$-valuesAn alternative way of doing the analysis above is to compute a $p$-value for an obtained experimental result.- Assuming that $H_0$ is true and the corresponding distribution of experimental results is known we accept some significance level $\alpha$. Usually $\alpha=0.05$.- Given an experimental result $x$ we compute the probability $P$ to see a value at least as extreme as $x$. - The extreme means that it deviates from the center of the distribution more then $x$. - For a one-sided test (e.g., weather a coin is head biased) this probability is $p$-value, $p=P$.- For two sides test (weather a coin is unfair) $p$-value is doubled this probability, $p = 2 P$.- If $p<\alpha$ the $H_0$ is rejected.To illustrate it, let us consider again a test if a coin is unfair. The null hypothesis $H_0$ is that the coin is fair. The experiment is tossing the coin $N$ times and the result is  a number of heads $x$.For this experiment the distribution of results $x$ is known: For large $N$ this is normal distribution.Let the significance level be $\alpha=0.05$.We conduct the experiment, compute $p$-value and reject $H_0$ if $p<\alpha$.<jupyter_code>import numpy as np
from scipy.special import erf

def get_mu_sig(N, P):
    """
    mu and sig corresponding to a binomial distribution"""
    mu = P * N
    sig = np.sqrt(N * P * (1 - P))
    return mu, sig

def norm_cdf(x, mu, sig):
    """Normal cumulative distribution function"""
    return 0.5 * (1 + erf((x-mu)/(sig*np.sqrt(2))))

def get_pvalue_twosided(x, mu0, sig0):
    if x <= mu0:
        return 2 * norm_cdf(x, mu0, sig0)
    else:
        return 2 * (1 - norm_cdf(x, mu0, sig0))

N = 1000
P = 0.5
alpha = 0.05

# mean and standard deviation for a fair coin
mu0, sig0 = get_mu_sig(N, P)

# Let the result of the experiment be
x = 530
p = get_pvalue_twosided(x, mu0, sig0)
print(f'x={x}, p-value={p}, alpha={alpha}, {"reject H0" if p < alpha else "H0 is true"}')

# Let the result of the experiment be
x = 531
p = get_pvalue_twosided(x, mu0, sig0)
print(f'x={x}, p-value={p}, alpha={alpha}, {"reject H0" if p < alpha else "H0 is true"}')<jupyter_output><empty_output><jupyter_text>In fact we do the same think as above. But now we do not compute the bounds `X_lo` and `X_hi` in advance.Instead we compute $p$-value for the experiment result $x$ and test if $H_0$ can be rejected comparing the $p$-value with the significance $\alpha$.Since the level of the significance is commonly known often just a $p$-value of the discovered feature is reported. Once again: $p$-value is the probability to see $x$ and more extreme values provided that $H_0$ is true. If this probability is too small we can conclude that this could not happen at all so that $H_0$ is not true and have to be rejected.  Confidence IntervalAnother way of treating the result of the stochastic experiment is computing a confidence interval. Assume again that we have a coin and want to know if it is unfair. We have performed a series of $N=1000$ tossing and get $x$ heads. From this result we can estimate the probability of this coin as $$\hat P=x/N$$If $N\to \infty$ our estimate $\hat P$ approaches true probability of the coin $P$.But how $\hat P$ deviates from $P$ for a finite $N$? The answer is given by a confidential interval.The confidence interval is the interval around the experimental value $\hat P$ such that the true value $P$ belongs to it with  the probability $1-\alpha$, where $\alpha$ is as usual the significance. Often $1-\alpha=0.95$.We need to estimate the distance between true value $P$ and its  approximation $\hat P$. When the series length $N$ is large, say $N>100$, it can be as follows: First we need to know the width of the distribution of $\hat P$. Let us remember that if $z$ is the standard normal distribution ($\mu=0$, $\sigma=1$) then $x$ is distributed with $\mu$ and $\sigma$ if$$x = \sigma z + \mu$$Divide it by $N$:$$x/N = (\sigma z + \mu) / N$$Since $\hat P = x/N$ we have:$$\hat P = \left(\frac{\sigma}{N}\right) z + \left(\frac{\mu}{N}\right)$$So, if the standard deviation of $x$ is $\sigma$ then the standard deviation for $\hat P$ is $\sigma/N$.Also if $x$ has a mean $\mu$ then the mean for $\hat P$ is $\mu/N$.Let us now remember that if $x$ is the number of heads of a coin with the probability $P$ the standard deviation is$$\sigma = \sqrt{N P (1-P)}, \mu=P N$$But we do not know $P$ and use $\hat P$ instead (do not forget that we know $\hat P$ since this is the value that we have obtained in the experiment).So instead we take $$\sigma' = \sqrt{N \hat P (1- \hat P)}, \mu=\hat P N$$ Finally we obtain an estimate for the unknown standard deviation for $\hat P$:$$\hat \sigma = \frac{\sigma'}{N}=\sqrt{\hat P (1- \hat P)/N}, \hat \mu = \hat P$$This is where we employ the assumption that $N$ is large. For small $N$ we cannot take $\hat \sigma$ as an approximate estimate of an known distribution width and must use another approach.Now the confidential interval is found as the $(1-\alpha)$ probability range, e.g., 95% range, for the normal distribution with $\mu=\hat \mu$ and $\sigma = \hat\sigma$.For $\alpha=0.05$ this is 95% range that can be found as$$[\hat \mu - 1.96 \hat \sigma, \hat \mu + 1.96 \hat \sigma]$$or$$\left[\hat P - 1.96 \sqrt{\hat P (1- \hat P)/N},\hat P + 1.96 \sqrt{\hat P (1- \hat P)/N}\right]$$Often the deviations $\pm 1.96 \hat \sigma$ are treated as an error range of the experimental result. We can say that true probability $P$ is obtained as$$P = \hat P \pm 1.96 \hat \sigma$$ Lets us do some checks. Below is the function that computes the confidential interval according to the formula above<jupyter_code>import numpy as np

def conf_interv_95(N, P):
    """Compute confidential interval.
    """
    return 1.96 * np.sqrt((P * (1 - P)) / N)<jupyter_output><empty_output><jupyter_text>Assume some results and compute the confidential interval for it:<jupyter_code># This is our result: number of heads after N tossing
x = 525

N = 1000

P_hat = x / N
ci = conf_interv_95(N, P_hat)
print(f'Estimated probability is P_hat={P_hat}±{ci:6.4f}')
print(f'Confidential interval is [{(P_hat - ci):6.4f}-{(P_hat + ci):6.4f}]')<jupyter_output><empty_output><jupyter_text>The meaning of all this computations is as follows:Thus we have conduct an experiment that includes $N=1000$ tossing of a coin whose probability $P$ is unknown. The result is $x=525$ of heads and thus the probability is estimated as $\hat P = 0.525$.The error of this measurements is $\pm 0.0310$.The confidential interval, i.e., a range where the true value $P$ falls with the 95% probability is $[0.4940, 0.5560]$.In particular it means that the coin still can be fair, because $P=0.5$ belongs to this interval.<jupyter_code># The coin stil can be fair
x = 530

N = 1000

P_hat = x / N
ci = conf_interv_95(N, P_hat)
print(f'Estimated probability is P_hat={P_hat}±{ci:6.4f}')
print(f'Confidential interval is [{(P_hat - ci):6.4f}-{(P_hat + ci):6.4f}]')
# Now we conclude that the coin is unfair
x = 531

N = 1000

P_hat = x / N
ci = conf_interv_95(N, P_hat)
print(f'Estimated probability is P_hat={P_hat}±{ci:6.4f}')
print(f'Confidential interval is [{(P_hat - ci):6.4f}-{(P_hat + ci):6.4f}]')
# If this value is obtained, we again conclude that the coin is unfair
x = 469

N = 1000

P_hat = x / N
ci = conf_interv_95(N, P_hat)
print(f'Estimated probability is P_hat={P_hat}±{ci:6.4f}')
print(f'Confidential interval is [{(P_hat - ci):6.4f}-{(P_hat + ci):6.4f}]')<jupyter_output><empty_output><jupyter_text>Two value above for which the coin can be treated as unfair are 469 and 531. Observe that this is exactly the values that we have obtained above performing the hypothesis testing.  A/B testIn electronic commerce a conversion rate is a portion of site visitors who did some activity. For example, during a week there was $N$ visitors and $n$ of them have bought something. Or among $N$ visitors $n$ have clicked an advertisement.Thus the conversion rate is$$c = n / N$$Assume that we have a web site and its design is changed. The question is if the new design brings more customers.Let us denote the old design as A and the new one as B. The A/B test is about testing of user preferences. To perform this test we start showing two versions of the website at random: some sees A and others obtain B.After a while we count $N_A$ visits of the website A and $N_B$ visitors of B.Among those who get A there are $n_A$ visitors who did something useful, and for B-visitors this number is $n_B$.The conversion rates for two version of the website are$$c_A = n_A / N_A, c_B = n_B / N_B$$These conversion rates are random variables. In fact we have only single values for them, but if we repeat the measurements many times we will obtain different $c_A$ and $c_B$. If $N_A$ and $N_B$ are large it is reasonable to treat them as having the normal distributions.Also we assume that $c_A$ and $c_B$ are independent. This is indeed the case if each visitor has access to the only one version of the website.As we already discussed above mean and standard deviations for $c_A$ and $c_B$ can be estimated as$$\mu_A = c_A, \sigma_A = \sqrt{c_A(1-c_A)/N_A}$$$$\mu_B = c_B, \sigma_B = \sqrt{c_B(1-c_B)/N_B}$$The null hypothesis for our analysis is that $\mu_A=\mu_B$, i.e., for infinitely many visitors $c_B = c_A$ and there is no difference between two versions of the website.Consider a variable $$c_{AB} = c_B - c_A$$Again, we have only one value of it. But repeating the measurements many times we will obtain different values. It can be derived that $c_{AB}$ is also normally distributed with the mean and the standard deviation$$\mu_{AB} = \mu_B - \mu_A, \sigma_{AB} = \sqrt{\sigma_A^2 + \sigma_B^2}$$According to our null hypothesis, $\mu_{AB}=0$.To test it we compute $p$-value for $c_{AB}$.Let us remember that it means that provided that the null hypothesis is true we need to computethe probability to see a value at least as extreme as $c_{AB}$.If $c_{AB}>0$ we compute it as follows:$$p = 2(1 - \phi(c_{AB}, \mu_{AB}=0, \sigma_{AB}))$$and for $c_{AB}<0$ we do$$p = 2\phi(c_{AB}, \mu_{AB}=0, \sigma_{AB})$$where $\phi(\cdot)$ is CDF (cumulative distribution function) for a normal distribution that corresponds to our null hypothesis.Notice that we multiply the probability by 2. This is because we perform a two sided test: Deviation to both sided are meaningful. If the $p$-value is below some accepted significance $\alpha$ level we conclude that the new website results in the change of the conversion ratio, i.e., $H_0$ is rejected.There two cases: If $c_{AB}>0$ the new website works better then the old one, and for $c_{AB}<0$ we conclude that the new website make things worse.Let us consider a particular example.<jupyter_code>import numpy as np
from scipy.special import erf

# The website A
Na = 1000
na = 200

# The website B
Nb = 1200
nb = 290

# Convertion ratios
ca = na / Na
cb = nb / Nb

cab = cb - ca

def sigma_test(N, c):
    """Standard deviation for c"""
    return np.sqrt((c * (1 - c)) / N)

def norm_cdf(x, mu, sig):
    """Normal cumulative distribution function"""
    return 0.5 * (1 + erf((x-mu)/(sig*np.sqrt(2))))

sig_a = sigma_test(Na, ca)
sig_b = sigma_test(Nb, cb)
sig_ab = np.sqrt(sig_a**2 + sig_b**2)
mu_ab = 0  # this accirding to our null hypothesis

p = 2*norm_cdf(cab, mu=mu_ab, sig=sig_ab) if cab < 0 else 2*(1 - norm_cdf(cab, mu=mu_ab, sig=sig_ab))

print(f"ca={ca:6.4f}, cb={cb:6.4f}, cab={cab:6.4f}, p-value={p:6.4f}")<jupyter_output><empty_output><jupyter_text>In this example $c_b$ is higher then $c_a$, but the difference is not very large.Do not forget that the values are volatile, so that we can not trust just to a mere positive difference.But the $p$-value is $\approx 0.02$. It means that we have to accept the new design since the significance of the hypothesis that it is better is very high. Even higher than usual level $0.05$.Let us remember: this $p$ value means that the probability that the new design has produced higher conversion ration by chance is $0.02$, that is very small.  Bayesian inferenceLet us remember Bayes' theorem.First we know that the probability to meet a person with a certain disease $D$ is $P(D)$. This is the prior probability.Then we get some test system that mines new data. It can be said that it provides an evidence for $D$.The test system operates as follows: It gives a positive result $T$ for a really diseased person with the probability$P(T|D)$ and erroneously detects the disease for a healthy person with the probability $P(T|\bar D)$.Using the new data produced by this test system we can refine the prior probability using Bayes' theorem:The probability to meet a diseased person provided that he or she already have a positive test is$$P(D|T) = \frac{P(T|D) P(D)}{P(T|D) P(D) + P(T|\bar D) P(\bar D)}$$This probability is the posterior probability after taking the evidence $T$ into account$P(T|D)$ is called a likelihood function. This is the probability of the evidence $T$ provided that $D$ is true.Thus the Bayes' theorem refines somehow known prior probability using a new evidence.  Bayes' theorem can be used in statistical inference exactly as shown above.First we have some prior estimates of a probability distribution. This can be even intuitive guesses or expert assessments.Then we find a new information and refine the prior probabilities via Bayes' theorem.Using Bayes' theorem for continuous variables is rather complicated mathematical problem: Given a prior probability distribution function we have to find a posterior distribution. But if a prior distribution function is prior conjugate to a likelihood distribution the posterior distribution function is the same as the prior one but with another parameters.Many prior conjugate pairs of distributions are known.The experimental results, i.e., the evidences considered in the examples above (number of heads in tossing of a coin and number of consumers among all website visitors) have binomial distribution. The conjugate pair for it is so called beta distribution.Before proceeding we will consider what is it.Beta distribution is$$\rho(x) = \frac{ x^{a-1} (1-x)^{b-1} }{B(a,b)}$$where $B(a,b)$ is a beta function:$$B(a,b)=\int_0^1 x^{a-1}(1-x)^{b-1} dx$$In fact we will not use the former formula to compute $B(a,b). The module `scipy.special` provides a function `beta` for it.<jupyter_code>from scipy.special import beta
import numpy as np

def beta_pdf(x, a, b):
    if x <= 0 or x >= 1:
        return 0
    return x**(a - 1) * (1 - x)**(b - 1) / beta(a, b)<jupyter_output><empty_output><jupyter_text>First of all notice that this distribution is located in the unit interval and is zero everywhere outside it.There are two parameters, $a$ and $b$. For $a=b=1$ it is just the uniform distribution. When $a=b>1$ the curve is symmetric with respect to $1/2$ and has a hump-like shape. The higher values the more pointy and narrow is the curve. At $ab$ all probabilities are near 1.The mean value for this distribution is$$\mu=\frac{\alpha}{\alpha+\beta}$$<jupyter_code>import matplotlib.pyplot as plt

xx = np.linspace(0, 1, 100)
y_1_1 = [beta_pdf(x, 1, 1) for x in xx]
y_10_10 = [beta_pdf(x, 10, 10) for x in xx]
y_30_30 = [beta_pdf(x, 30, 30) for x in xx]
y_5_10 = [beta_pdf(x, 5, 10) for x in xx]
y_20_5 = [beta_pdf(x, 20, 5) for x in xx]

fig, ax = plt.subplots()
ax.plot(xx, y_1_1, label='a=1, b=1')
ax.plot(xx, y_10_10, label='a=10, b=10')
ax.plot(xx, y_30_30, label='a=30, b=30')
ax.plot(xx, y_5_10, label='a=5, b=10')
ax.plot(xx, y_20_5, label='a=20, b=5')
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>As already mentioned above the beta distribution is conjugate prior to the binomial distribution. It means that if a priory distribution is beta we can refine it using an experiment describing by the binomial distribution.It works as follows. First we choose $a$ and $b$ to fit the priory probability $P$ of a considered Bernoulli variable, e.g., a coin.$$P = \frac{a}{a+b}$$If there are no guess about $P$ at all we set $a=b=1$ that corresponds to the uniform distribution. If there are reasons to expect that $P\approx 0.5$ we set $a=b>1$. Larger values of $a$ and $b$ correspond to higher confidence in this guess since for larger $a$ and $b$ the curve becomes more concentrated near $0.5$.For other guesses of $P$ we can also find $a$ and $b$ in a similar way: the mean value of the distribution corresponds to the guess and its width reflects our confidence in it.Then we perform an experiment, say toss a coin. Let as usual there $N$ tossing and $x$ heads. According to the Bayes' theorem the posterior distribution, i.e., the one refined after new obtaining the new information will also be beta with $$a'=a+x, b'=b+(N-x)$$The mathematical details are too complicated to be mentioned here.Thus a refined estimate of $P$ is$$P' = \frac{a+x}{a+b+N}$$Below is an example of computations using the above formula.<jupyter_code>def bayes_infer(a, b, N, x):
    """Takes initial a and be and computes the refined a1 and b1.
    Returns priori and posterir P.
    """
    P = a / (a + b)
    a1 = a + x
    b1 = b + (N - x)
    P1 = a1 / (a1 + b1)
    return P, P1

a = 1
b = 1
N = 1000
x = 530

P, P1 = bayes_infer(a, b, N, x)
print(f"prior P={P:6.4f}, posterior P1={P1:6.4f}")<jupyter_output><empty_output><jupyter_text>Example of using of Bayesian inferenceConsider a more realistic example. Assume that there is a website and previously we have computed its conversion ratio as$$c = 0.11$$Now we have changed its design and want to refine the conversion ratio. The initial guess is obviously our previous conversion rate. We need to find such $a$ and $b$ that$$c = \frac{a}{a+b}$$There is one equation for two unknown variables, $a$ and $b$. Solve it for $b$:$$b = a \frac{1-c}{c}$$We have one free variable $a$: choosing it we compute the corresponding $b$.Let us see how the distribution depends on $a$:<jupyter_code>import matplotlib.pyplot as plt

c = 0.11

xx = np.linspace(0, 0.3, 100)

fig, ax = plt.subplots()

for a in [1,2,5]:  # possible values of a
    b = a * (1-c)/c  # corresponding b
    yy = [beta_pdf(x, a, b) for x in xx]
    ax.plot(xx, yy, label=f'a={a:6.4f}, b={b:6.4f}')
    
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>The width of the curve reflects to our expectations.The higher $a$ results the narrower curve. It corresponds to the lower uncertainty in our expectations: we expect that the new design will have rather similar performance as the old one and $c$ will be changed slightly.But if we expect that the new design will result in the essential corrections, we have to keep $a$ as low as possible. Let us be optimistic  and expect a serious correction. At $a=1$ the curve looks inappropriately, so we set$$a = 2, b = 16.1818$$Now we start the test: count the total number of visitors $N$ and those who did something useful, e.g., bought something. The former will be denoted as $x$.Let us assume that after some time we have$$N = 100, x = 17$$Using the Python function defined above we refine the conversion ratio as follows<jupyter_code>c = 0.11
a = 2
b = a * (1 - c) / c

N = 100
x = 17

_, c1 = bayes_infer(a, b, N, x)

print(f"prior c={c:6.4f}, posterior c1={c1:6.4f}, plain mean x/N={x/N:6.4f}")<jupyter_output><empty_output><jupyter_text>The updated conversion ratio is$$c'=0.1608$$It indicates that the improvements have occurred indeed.Notice that the new value $c'$ differs from the experimental mean value $$x/N = 0.17$$Here is the plots of the prior and posterior distributions. Observe that the posterior one is narrower: It reflects that fact that after obtaining a new information the uncertainty is decreased.<jupyter_code>import matplotlib.pyplot as plt

c = 0.11
a = 2
b = a * (1 - c) / c

N = 100
x = 17

a1 = a + x
b1 = b + (N - x)

xx = np.linspace(0, 1, 500)

fig, ax = plt.subplots()
y0 = [beta_pdf(x, a, b) for x in xx]
y1 = [beta_pdf(x, a1, b1) for x in xx]
    
ax.plot(xx, y0, label=f'a={a:6.4f}, b={b:6.4f}')
ax.plot(xx, y1, label=f'a1={a1:6.4f}, b1={b1:6.4f}');
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>This analysis makes sense only if we have a small amount of data.Above there were $N=100$ visitors.Let us compute $c'$ if $N=1000$:<jupyter_code>c = 0.11
a = 2
b = a * (1 - c) / c

# incresed by a factor of 10
N = 1000
x = 170

_, c1 = bayes_infer(a, b, N, x)

print(f"prior c={c:6.4f}, posterior c1={c1:6.4f}, plain mean x/N={x/N:6.4f}")<jupyter_output><empty_output><jupyter_text>We see that now $c'$ is very close to the experimental mean $x/N$.Below are the corresponding plots:<jupyter_code>import matplotlib.pyplot as plt

c = 0.11
a = 2
b = a * (1 - c) / c

# incresed by a factor of 10
N = 1000
x = 170

a1 = a + x
b1 = b + (N - x)

xx = np.linspace(0, 1, 500)

fig, ax = plt.subplots()
y0 = [beta_pdf(x, a, b) for x in xx]
y1 = [beta_pdf(x, a1, b1) for x in xx]
    
ax.plot(xx, y0, label=f'a={a:6.4f}, b={b:6.4f}')
ax.plot(xx, y1, label=f'a1={a1:6.4f}, b1={b1:6.4f}');
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>The distribution is very pointy. It tells us that we must be pretty sure that the new conversion ration is $x/N=0.17$.Thus we conclude that when a lot of experimental data is available the Bayesian inference is not so useful since its result mere coincides with the experimental mean.   Chi-square distributionLet us remember that many natural processes can be modeled as stochastic ones with a normal distribution. This due to the central limit theorem that says that if many independent factors influences an observable then their impacts add up so that the observable demonstrates a random behavior with a normal distribution.The PDF (probability density function) for the normal distribution is as follows $$\rho(x)=\frac{1}{\sigma\sqrt{2\pi}}\mathrm{exp}\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$When $\mu=0$ and $\sigma=1$ the distribution is called the standard normal distribution.The plot of this distribution is as follows:<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

def norm_pdf(x, mu, sig):
    """Normal probability density function"""
    return np.exp(-(x-mu)**2 / (2*sig*sig)) / (sig * np.sqrt(2*np.pi))

mu, sig = 0, 1

x = np.linspace(-5, 5, 100)
y = [norm_pdf(xi, mu, sig) for xi in x]

fig, ax = plt.subplots()

ax.plot(x, y, label=f"mu={mu}, sig={sig}")

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho$')
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>Sum and difference of two random normal variables have also the normal distribution. Now let us remember what is the distance. Assume that we have two vectors. $$\vec x = (x_1, x_2, x_3, x_4, x_5)$$$$\vec y = (y_1, y_2, y_3, y_4, y_5)$$In fact these are just two sets of numbers.If we want to compare how different they are we need a distance between them. Most often the Euclidean distance is used:$$d = \sqrt{\sum_{i=1}^{k} (x_i - y_i)^2}$$where $k$ is the number of vector components. For two vectors above $k=5$.Often comparing two vectors the square root is omitted. The meaning of the result is not changed, but computations are simplified.$$d^2 = \sum_{i=1}^{k} (x_i - y_i)^2$$ Now assume that we have two random vectors with $k$ elements and the elements are random normal variables.$$\vec x = (x_1, x_2, x_3, x_4, x_5)$$$$\vec y = (y_1, y_2, y_3, y_4, y_5)$$Often we need to answer the question: Do elements of these two vectors belong to the same probability distribution or not?The answer is given with the help of $\chi^2$-test. Consider $k$ independent random numbers that have standard normal distribution, i.e., $\mu=0$, $\sigma=1$:$$z_1, z_2, z_3, z_4, z_5$$The distribution of the sum of their squares$$d^2 = \sum_{i=1}^k z_i^2$$is called $\chi^2$-distribution. This distribution has one parameter $k$ that is a number of $z$'s. This $k$ is called a number of degrees of freedom.PDF for $\chi^2$-distributions is given by a formula:$$\rho(x) = \frac{(1/2)^{\frac{k}{2}}}{\Gamma\left(\frac{k}{2}\right)} x^{\frac{k}{2} - 1} e^{-\frac{x}{2}}$$where $\Gamma(\cdot)$ is gamma function. We will compute if using function `gamma` from `scipy.special`.Examples of the $\chi^2$ PDF for different $k$ are shown below.<jupyter_code>from scipy.special import gamma
import numpy as np

def chi2_pdf(x, k):
    """Chi-square PDF"""
    return (1/2)**(k/2) * x**(k/2-1) * np.exp(-x/2) / gamma(k/2)
import matplotlib.pyplot as plt

xx = np.linspace(1e-5, 10, 100)
y_1 = [chi2_pdf(x, 1) for x in xx]
y_2 = [chi2_pdf(x, 2) for x in xx]
y_3 = [chi2_pdf(x, 3) for x in xx]
y_4 = [chi2_pdf(x, 4) for x in xx]
y_5 = [chi2_pdf(x, 5) for x in xx]

fig, ax = plt.subplots()
ax.plot(xx, y_1, label='k=1')
ax.plot(xx, y_2, label='k=2')
ax.plot(xx, y_3, label='k=3')
ax.plot(xx, y_4, label='k=4')
ax.plot(xx, y_5, label='k=5')
ax.set_ylim([0, 0.6])
ax.legend()
ax.grid()<jupyter_output><empty_output><jupyter_text>So, we want to compare two random vectors with $k$ elements and the elements are random normal variables. $$\vec x = (x_1, x_2, x_3, x_4, x_5)$$$$\vec y = (y_1, y_2, y_3, y_4, y_5)$$Though this test can be performed for an arbitrary sets of number, typically the compared vectors are frequencies of certain values of some stochastic variable. Thus the vector elements are assumed to be positive.In principle to compare $\vec x$ and $\vec y$ we could compute their plan distance and analyse its distribution. But it would be rather complicated problems since the original data can have an arbitrary magnitudes. Instead we compute the rescaled squared distance as$$d^2 = \sum_{i=1}^{k} \frac{(x_i - y_i)^2}{x_i}$$Then we select a significance $\alpha$, say $\alpha=0.05$, and check if the computed $d^2$ falls within $(1-\alpha)$ probability area of $\chi^2$-distribution. If yes the considered vectors are sampled from the same distributions.Typically this test is performed to compare the theoretically expected frequencies with the observed ones.Notice that previously we defined a frequency as number of a particular observed values divided by the total number of observations. They were number within the range between 0 and 1. More precisely it can be called a relative frequency.Performing $\chi^2$ test we call the number of observed values as the frequency. The frequency is an integer number. The sum of all frequencies equals to the total number of observations. Thus to perform $\chi^2$-test we compute$$d^2 = \sum_{i=1}^{k} \frac{(E_i - O_i)^2}{E_i}$$where $E_i$ are the expected frequencies and $O_i$ are the observed frequencies.Then this value is compared with the $\chi^2$-distribution as mentioned above.  Chi-square test of an unfair coinLet us again perform a test of a coin. The idea is that we conduct a series of tossing many times. Now $N$ is small, say $N=5$ and we repeat it $M=250$ times.The result of each series is a binomial variable $x$ within the range from 0 to $N$. Thus we have a series of $M$ generation of a binomial variable $x$. We count the frequencies of of each value of $x$.Also we take the binomial distribution and compute theoretically expected frequencies.Let us first find the theoretical frequencies corresponding to the $H_0$ hypothesis.<jupyter_code>import numpy as np
from scipy.special import binom

def binom_distr(n, k, p):
    """Binomial distribution"""
    return binom(n, k) * p**k * (1-p)**(n-k)

N = 5
M = 250

P_h0 = 0.5  # This P corresponds to the H0 hypothesis

freq_expected = [M * binom_distr(N, k, P_h0) for k in range(N+1)]

# expected frequencies
print(freq_expected)

# their sum equals to the total number of experiments M
print(sum(freq_expected))<jupyter_output><empty_output><jupyter_text>Here we model an experiment: repeat a series of $N$ tossing $M$ times and count observed frequencies.<jupyter_code>import numpy as np

P_actual = 0.55 # Actual probability of the coin

freq_observed = [0] * (N + 1) 
rng = np.random.default_rng()
for _ in range(M):  # repeat series of tossing M times
    # one series: toss a coin N times and freq_observed heads
    x  = sum([1 if x < P_actual else 0 for x in rng.random(size=N)]) 
    # freq_observed[0] - how many times zero heads appearer in a series of N tossing
    # freq_observed[1] - how many times one head appeared after N tossing 
    # freq_observed[2] - two head in N tossing 
    # and so on
    freq_observed[x] += 1  # freq_observed appeard number of heads

# computed freq_observeders
print(freq_observed)

# sum of all freq_observed must be M
print(sum(freq_observed))<jupyter_output><empty_output><jupyter_text>Now we perform the $\chi^2$-test. While the common idea of this test is not so complicated, the actual computations includes some subtle steps and we do not write the routine ourself. (For example frequencies with less then 5 observations must be merged with other frequencies.)The function `chisquare` from `scipy.stats` is used instead.<jupyter_code>from scipy.stats import chisquare

chi = chisquare(f_obs=freq_observed, f_exp=freq_expected)

# Print the result return by chisquare
print(chi)

alpha = 0.05

if chi.pvalue < alpha:
    print("H0 is rejected, the coin is unfair")
else:
    print("H0 is not rehected, the coin is fair")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11094475220841443, 'fraction_numerical': 0.024854020062883665, 'mean_word_length': 4.416421692853523, 'pattern_counts': {'":': 0, '<': 174, '<?xml version=': 0, '>': 160, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '37', 'dir': 'jupyter-structured-clean-dedup', 'id': '451386', 'license': 'MIT', 'n_tokens_mistral': 18131, 'n_tokens_neox': 17137, 'n_words': 7475, 'path': '06_Hypothesis_and_inference.ipynb', 'repo_name': 'kupav/data-sc-intro', 'size': '53450'}