<jupyter_start><jupyter_text>Day 2 Intro to Probability and Statistics  Probability: Expectation, Standard Deviation  Expected Value$(\Omega, \mathcal{F}, P)$ - [probability space](https://en.wikipedia.org/wiki/Probability_space)Let $\xi$ - some [random variable](https://en.wikipedia.org/wiki/Random_variable): $\xi: \Omega \rightarrow X$ (Usually $X \subseteq \mathbf{R}$). Let $X = \{x_1,\dots x_k\}$and $A_i = \{w\in\Omega: \xi(w) = x_i \}$Then, you can represent $\xi$ as follows:$\xi(w) = \sum_{i=1}^k{x_i * I(A_i)}$where $I$ - [indicator function](https://en.wikipedia.org/wiki/Indicator_function):$I(A_i) = \left\{\begin{array}{ll}1,\ w \in A_i \\ 0,\ w \notin A_i \\\end{array}\right.$Then [Expected Value](https://en.wikipedia.org/wiki/Expected_value) $E$ is:$E(\xi) = \sum_{i=1}^k{x_i * P_\xi(A_i)}$where $P_\xi$ is [probability distribution](https://en.wikipedia.org/wiki/Probability_distribution)  Expected value: Properties1. If $\xi \geq 0 \Rightarrow E(\xi) \geq$2. $E(a\xi + b\eta) = aE(\xi) + bE(\eta)$3. If $\xi \geq \eta \Rightarrow E(\xi) \geq E(\eta)$4. $|E(\xi)| \leq E(|\xi|)$5. $\xi \bot \eta \Rightarrow E(\xi\eta) = E(\xi)E(\eta)$6. $(E|\xi\eta|)^2 \leq E(\xi^2)E(\eta^2)$ - [Cauchy–Schwarz inequality](https://en.wikipedia.org/wiki/Cauchy–Schwarz_inequality)7. $\xi = I(A) \Rightarrow E(\xi) = P(A)$  Variance[Variance](https://en.wikipedia.org/wiki/Variance) of random variable $\xi$ is the expectation of the squared deviation of $\xi$:$D(\xi) = E(\xi - E(\xi))^2$[Standard deviation](https://en.wikipedia.org/wiki/Standard_deviation): $\sigma = +\sqrt{D}$  Variance: Properties1. $D(\xi) = E(\xi - E(\xi))^2 = E(\xi^2 - 2\xi E(\xi) + (E(\xi))^2) = E(\xi^2) - 2(E(\xi))^2 + (E(\xi))^2 = E(\xi^2) - (E(\xi))^2$2. $D(a) = 0$3. $D(\xi) = 0 \leftrightarrow\exists a: P(\xi = a) = 1$ 4. $D(a\xi) = a^2D(\xi)$5. $D(\xi + b) = D(\xi)$6. $\xi\bot\eta\Rightarrow D(\xi+\eta)=D(\xi)+D(\eta)$*Question:* $D(\xi + \eta) = ?$ --- Okay, let see some examples and try to empirically confirm some properties<jupyter_code>import scipy.stats as sps
import numpy as np<jupyter_output><empty_output><jupyter_text>Example: [Bernoulli distribution](https://en.wikipedia.org/wiki/Bernoulli_distribution) $E(X) = \sum_i{x_i * p(X = x_i)}$$X = \left\{\begin{array}{ll}1,\ p \\ 0,\ 1 - p \\\end{array}\right.$$E(X) = 1 * P(X = 1) + 0 * P(X = 0) = P(X = 1)$---$$D(X) = E(X^2) - (E(X))^2 = P(X=1) - P(X=1)^2 = P(X=1)*(1 - P(X=1)) = P(X=1)P(X=0)$$<jupyter_code>p = 0.5
size = 100000

x = sps.bernoulli.rvs(p, size=size)
# Very bad way
# Only for studying purposes
def get_value_prob(value: int, x: np.ndarray) -> float:
        return np.sum(x == value) / len(x)


def calc_expectation(x: np.ndarray) -> float:
    E = 0.
    for value in set(x.tolist()):
        E += value * get_value_prob(value, x)
    return E


def calc_variance(x: np.ndarray) -> float:
    var = 0.
    E = calc_expectation(x)
    for value in set(x.tolist()):
        var += (value - E)**2 * get_value_prob(value, x)
    return var
print(f"Expected Value: {calc_expectation(x)}\nVariance: {calc_variance(x)}")<jupyter_output><empty_output><jupyter_text>Expected Value: Linear property + Multiplication of independent random variables<jupyter_code>a = 3
b = 4
p1 = 0.1
p2 = 0.3
size = 100000
x1 = sps.bernoulli.rvs(p1, size=size)
x2 = sps.bernoulli.rvs(p2, size=size)
x = a * x1 + b * x2

print(f"E(X) = {calc_expectation(x)}"\
     f"\naE(X1) + bE(X2) = {a*calc_expectation(x1) + b*calc_expectation(x2)}"\
     f"\na*p1 + b*p2 = {a*p1 + b*p2}"
)
print("-"*25)
print(f"E(X1*X2) = {calc_expectation(x1 * x2)}"\
     f"\nE(X1)E(X2) = {calc_expectation(x1)*calc_expectation(x2)}")<jupyter_output><empty_output><jupyter_text>Variance: properties<jupyter_code>print(f"D(aX1) = {calc_variance(a * x1)}\n"\
      f"a**2D(X1) = {a**2 * calc_variance(x1)}")
print("-"*25)
print(f"D(X1) = {calc_variance(x1)}\nD(X1 + b) = {calc_variance(x1 + b)}")
print("-"*25)
print(f"D(X1+X2) = {calc_variance(x1+x2)}"\
      f"\nD(X1)+D(X2) = {calc_variance(x1) + calc_variance(x2)}")
print("-"*25)
print(f"D(X1 + sin(X1)) = {calc_variance(x1 + np.sin(x1))}\n"\
      f"D(X1) + D(sin(X1)) = {calc_variance(x1) + calc_variance(np.sin(x1))}")<jupyter_output><empty_output><jupyter_text>Example: [Binomial distribution](https://en.wikipedia.org/wiki/Binomial_distribution)The binomial distribution with parameters n and p ($Bin(n,p)$) - the discrete probability distribution of the number of successes in a sequence of n independent experiments, each asking a yes–no question$EX = \text{(Linear property)} = \sum_i^n{EX_i}$where $X_i \sim Bern(p)$$EX = nEX_1 = np$---$$DX = (X_i \bot X_j\ \forall i, j \in [1\dots n], i \neq j) = \sum_i^n{DX_i} = nDX_1 = np(1 - p)$$<jupyter_code>n = 10
p = 0.5
x = sps.binom.rvs(n, p, size=size)

print(f"np = {n*p}; np(1 - p) = {n * p * (1 - p)}\n"\
      f"E(X) = {x.mean()}\n"\
      f"D(X) = {x.var()}"
)<jupyter_output><empty_output><jupyter_text>[Law of Large Numbers](https://en.wikipedia.org/wiki/Law_of_large_numbers)Let $\xi_1,\dots\xi_n$ - i.i.d; $E\xi_1 = a$$$\frac{\xi_1+\dots+\xi}{n} \rightarrow^{a.s.} a$$[Almost sure convergence](https://en.wikipedia.org/wiki/Convergence_of_random_variablesAlmost_sure_convergence) Let's try to visualize it<jupyter_code>import seaborn as sns
import matplotlib.pyplot as plt

sns.set(font_scale=1.6, palette='summer')
%matplotlib inline
size = 1000
samples = sps.bernoulli(p=0.5).rvs(size=size)
cum_means = samples.cumsum() / (np.arange(size) + 1)
plt.figure(figsize=(15, 5))
plt.plot(cum_means, lw=3)
plt.hlines(0.5, 0, size, alpha=0.3)
plt.xlabel('Number of random variables')
plt.ylabel('The value of the average')
plt.xlim((0, size));<jupyter_output><empty_output><jupyter_text>Single experiment is almost always not enough. Let's conduct 10 more experiments, and visualize each of them on different plots<jupyter_code>plt.figure(figsize=(20, 20))

for i in range(10):

    samples = sps.bernoulli(p=0.5).rvs(size=size)
    cum_means = samples.cumsum() / (np.arange(size) + 1)
    
    plt.subplot(5, 2, i+1)
    plt.plot(cum_means, lw=3)
    plt.hlines(0.5, 0, size, alpha=0.3)
    plt.xlabel('Number of random variables')
    plt.ylabel('The value of the average')
    plt.xlim((-5, size))
    
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>And now, even more experiments, on single plot<jupyter_code>size = 1000 
samples_count = 500

samples = sps.bernoulli(p=0.5).rvs(size=(samples_count, size))
cum_means = samples.cumsum(axis=1) / (np.arange(size) + 1)
plt.figure(figsize=(15, 7))

for i in range(samples_count):
    plt.plot(np.arange(size) + 1, cum_means[i], color='green', alpha=0.02)

plt.xlabel('Number of random variables')
plt.ylabel('The value of the average')
plt.xlim((0, size));<jupyter_output><empty_output><jupyter_text>[Central Limit Theorem](https://en.wikipedia.org/wiki/Central_limit_theorem) $X_1,\dots,X_n$ - i.i.d,    $$ \frac{\sum{X_i}}{n} \rightarrow^d \mathcal{N}\left(E(X), \frac{D(X)}{n}\right)$$    [Convergence in distribution](https://en.wikipedia.org/wiki/Convergence_of_random_variablesConvergence_in_distribution) Let's see an example.We take random variables from [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) and apply clt<jupyter_code>sample = sps.expon.rvs(size=1000)
fig, ax = plt.subplots(1, 1)

x = np.linspace(
    sps.expon.ppf(0.01),
    sps.expon.ppf(0.99), 
    100
)

ax.plot(x, sps.expon.pdf(x), 'r-', lw=5, alpha=0.5, label='Exponential pdf')
ax.hist(
    sample,
    bins=100,
    density=True,
    histtype='stepfilled',
    alpha=0.5,
    label='sample hist'
)

plt.xlabel('Count')
plt.ylabel('Frequency')
ax.legend(loc='best', frameon=False)
plt.show()
sample_sizes = [500, 1000, 5000]
sample_count = 1000
lambd = 2.

for size in sample_sizes:
    print(f"Sample size: {size}")
    samples = [
        sps.expon.rvs(size=size, scale=1./lambd).mean()
        for _ in range(sample_count)
    ]
    
    plt.hist(
        samples, bins=size//100, 
        density=True, histtype='stepfilled', 
        alpha=0.5, label='sample hist'
    )
    
    mu = lambd**(-1)
    sigma = (lambd**(-2)/size)**0.5
    x = np.linspace(min(samples),
                    max(samples), 100)
    
    plt.plot(x, sps.norm.pdf(x, loc=mu, scale=sigma),
             'r-', lw=5, alpha=0.6, label='norm pdf')
    
    plt.xlabel('Count')
    plt.ylabel('Frequency')
    plt.legend(loc='best', frameon=False)
    plt.show()<jupyter_output><empty_output><jupyter_text>[Statistical hypothesis testing](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing)  [Kolmogorov–Smirnov test](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test)$X_1\dots X_n$ - samples with unknown [cumulative distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function) $F$$H_0:F = F_0$$H_1:F = F_1$<jupyter_code>sps.kstest(sps.norm.rvs(size=100), cdf=sps.norm.cdf)
# sps.kstest?<jupyter_output><empty_output><jupyter_text>*KstestResult* has 2 fields`statistic` - criterion statistics`pvalue` - The p-value is the probability that a given result (or a more significant result) would occur under the null hypothesis.<jupyter_code>def reject_or_nor(pval: float, alpha: float = 0.05):
    print(f"P-value = {pval}")
    if pval < alpha:
        print(f"The null hypothesis is rejected at"
              f"the chosen level of significance {alpha}")
    else:
        print(f"The null hypothesis is NOT rejected at" 
              f"the chosen level of significance {alpha}")
# shouldn't be rejected:
pval = sps.kstest(sps.norm.rvs(size=100), cdf=sps.norm.cdf).pvalue
reject_or_nor(pval)

# should be rejected:
pval = sps.kstest(sps.norm(loc=10).rvs(size=100), cdf=sps.norm.cdf).pvalue
reject_or_nor(pval)


# should be rejected:
pval = sps.kstest(sps.uniform.rvs(size=100), cdf=sps.norm.cdf).pvalue
reject_or_nor(pval)<jupyter_output><empty_output><jupyter_text>[Shapiro-Wilk criterion](http://www.machinelearning.ru/wiki/index.php?title=Критерий_Шапиро-Уилка)$X_1\dots X_n$ - samples $H_0:$ samples are from [normal distribution](https://ru.wikipedia.org/wiki/Нормальное_распределение) with arbitrary parameters$H_1:$ samples aren't from normal distribution<jupyter_code># shouldn't be rejected:
pval = sps.shapiro(sps.norm.rvs(size=100)).pvalue
reject_or_nor(pval)

# shouldn't be rejected:
pval = sps.shapiro(sps.norm(loc=20, scale=100).rvs(size=100)).pvalue
reject_or_nor(pval)


# should be rejected:
pval = sps.shapiro(sps.uniform.rvs(size=100)).pvalue
reject_or_nor(pval)<jupyter_output><empty_output><jupyter_text>[Student's t-test (independent samples)](https://en.wikipedia.org/wiki/Student%27s_t-testIndependent_(unpaired)_samples)$X_1\dots X_n\sim\mathcal{N}(a_1, \sigma_1^2)$$Y_1\dots Y_m\sim\mathcal{N}(a_2, \sigma_2^2)$$X_i\bot Y_j$---$H_0: a_1 = a_2$$H_1: a_1 \neq a_2$<jupyter_code># should be rejected
sample_1 = sps.norm(loc=0).rvs(size=100)
sample_2 = sps.norm(loc=1).rvs(size=150)
pval = sps.ttest_ind(sample_1, sample_2).pvalue
reject_or_nor(pval)

# should be rejected, but sometimes....
sample_1 = sps.norm(loc=0).rvs(size=50)
sample_2 = sps.norm(loc=1, scale=7).rvs(size=100)
pval = sps.ttest_ind(sample_1, sample_2, equal_var=False).pvalue

reject_or_nor(pval)<jupyter_output><empty_output><jupyter_text>[Student's t-test (paired samples)](https://en.wikipedia.org/wiki/Student%27s_t-testPaired_samples)$X_1\dots X_n\sim\mathcal{N}(a_1, \sigma_1^2)$$Y_1\dots Y_m\sim\mathcal{N}(a_2, \sigma_2^2)$$X_i$ and $Y_j$ can be paired---$H_0: a_1 = a_2$$H_1: a_1 \neq a_2$<jupyter_code># shouldn't be rejected
sample_1 = sps.norm(loc=0).rvs(size=100)
sample_2 = sample_1 + sps.norm(loc=0, scale=0.5).rvs(size=100)
pval = sps.ttest_rel(sample_1, sample_2).pvalue
reject_or_nor(pval)

# should be rejected
sample_1 = sps.norm(loc=0).rvs(size=100)
sample_2 = sample_1 + sps.norm(loc=1, scale=0.5).rvs(size=100)
pval = sps.ttest_rel(sample_1, sample_2).pvalue
reject_or_nor(pval)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18881821225309675, 'fraction_numerical': 0.032725142283227315, 'mean_word_length': 4.679182509505703, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 60, 'https://': 21, 'lorem ipsum': 0, 'www.': 1, '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': '38059', 'license': 'Apache-2.0', 'n_tokens_mistral': 5185, 'n_tokens_neox': 4876, 'n_words': 880, 'path': '02_math_recap/day02_probability.ipynb', 'repo_name': 'Vitaly-lv/madmo-basic', 'size': '11948'}