<jupyter_start><jupyter_text>Exploratory Data Analysis using statistical analysis  (Prieskumná analýza s využitím štatistickej analýzy) Asumptions- the data must be randomly sampled from the population of interest - the data variables follow a normal distribution- a reasonably large sample size is used. - homogeneity of variance, i.e., the standard deviations of samples are approximately equal.<jupyter_code>import matplotlib.pyplot as plt
import pylab as py
import seaborn as sns
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.stats as sm_stats
import statsmodels.stats.api as sms
import scipy.stats as stats
import math
from numpy.random import seed
from numpy.random import rand
from numpy.random import randn
from numpy import mean
from numpy import var
from math import sqrt
from scipy.stats import norm<jupyter_output><empty_output><jupyter_text>Normal distribution = Gaussian distribution<jupyter_code>mu = 0
variance = 1
sigma = math.sqrt(variance)
x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
plt.plot(x, stats.norm.pdf(x, mu, sigma))
plt.show()<jupyter_output><empty_output><jupyter_text>Kurtosis<jupyter_code>x = np.arange(-3, 3, 0.001)
y1 = norm.pdf(x_axis, 0, 1) 
y2 = norm.pdf(x_axis, 0, 0.5) 
plt.plot(x, y1, x, y2)<jupyter_output><empty_output><jupyter_text>Skewness<jupyter_code>from scipy.stats import skewnorm 
x = np.linspace(0, 5, 100) 
     
# Varying positional arguments 
y1 = skewnorm .pdf(x, 4, 0.5, 1) 
y2 = skewnorm .pdf(x, 2, 1, 1) 
plt.plot(x, y1, x, y2, "r--")<jupyter_output><empty_output><jupyter_text>1. Normality Checks**A large fraction of the field of statistics is concerned with data that assumes that it was drawn from a Gaussian distribution (also called normal distribution).**  1.1 Histogram<jupyter_code>data_points = np.random.normal(0, 1, 100)

df = pd.DataFrame(data_points)
df.hist()
sns.histplot(df, bins=10)<jupyter_output><empty_output><jupyter_text>1.2 Distribution plotURL https://machinelearningknowledge.ai/seaborn-histogram-plot-using-histplot-tutorial-for-beginners/<jupyter_code># sns.distplot(df, bins=10)
sns.displot(df)<jupyter_output><empty_output><jupyter_text>1.3 Q-Q plot<jupyter_code>import statsmodels.api as sm

sm.qqplot(data_points, line='45')
py.show()<jupyter_output><empty_output><jupyter_text>**Q-Q plot interpretation**URL: https://stats.stackexchange.com/questions/101274/how-to-interpret-a-qq-plot  1.4 Shapiro-Wilk test of normality**The Shapiro-Wilk test tests the null hypothesis $H_0$ that the data was drawn from a normal distribution. The chance of rejecting the null hypothesis $H_0$ when it is true is close to 5% regardless of sample size.**URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html<jupyter_code>from scipy.stats import shapiro

np.random.seed(123)

x = stats.norm.rvs(loc=5, scale=3, size=100)
# x = stats.norm.rvs(loc=5, scale=3, size=1000)
# sns.distplot(x, bins=10)
# sns.displot(x, bins=10)
sns.histplot(x, bins=10)

shapiro_test = shapiro(x)
print(shapiro_test)

# interpret
alpha = 0.05
if shapiro_test.pvalue > alpha:
    print('Normal distribution (fail to reject H0)')
else:
    print('Another distributions (reject H0)')<jupyter_output><empty_output><jupyter_text>2. Statistical hypothesis testing**Do normality checks before hypothesis testing**- Histogram- Distribution plot- Q-Q plot- Shapiro-Wilkov Test  2.1 Levene test - Variance test**Hypothesis**- $H_0$: all input samples are from populations with equal variances - **Fail to Reject $H_0$**: Equal variances - **Reject $H_0$**: Another variancesURL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.levene.html<jupyter_code>from scipy.stats import levene
a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]

# calculate variances
[np.var(x, ddof=1) for x in [a, b, c]]
levene_test = levene(a, b, c)
print(levene_test)

# interpret
alpha = 0.05
if levene_test.pvalue > alpha:
    print('Equal variances (fail to reject H0)')
else:
    print('Another variances (reject H0)')<jupyter_output><empty_output><jupyter_text>**The small p-value < 0.05 suggests that the populations do not have equal variances.**  2.2 Student’s T-test $t = \frac{\mu_1 - \mu_2}{s_p \sqrt{\frac{1}{n_1} - \frac{1}{n_2}}}$ where $s_p = \sqrt{\frac{(n_1 - 1) s^2_{X_1} + (n_2 - 1) s^2_{X_2}}{n_1 + n_2  - 2}}$**Hypothesis**- $H_0$: the means of two populations are equal - **Fail to Reject $H_0$**: No difference between the sample means - **Reject $H_0$**: Some difference between the sample meansURL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html<jupyter_code>from scipy.stats import ttest_ind

seed(1)
# data1 = 5 * randn(100) + 50
# data2 = 5 * randn(100) + 51

data1 = stats.norm.rvs(loc=5, scale=10, size=500)
data2 = stats.norm.rvs(loc=5, scale=10, size=500)

# sns.distplot(data1, bins=10)
# sns.distplot(data2, bins=10)

sns.histplot(data1, bins=10, color='blue')
sns.histplot(data2, bins=10, color='green')

# compare samples
stat, p = ttest_ind(data1, data2)
print('Statistics=%.3f, p=%.3f' % (stat, p))

# interpret
alpha = 0.05
if p > alpha:
    print('Same distributions (fail to reject H0)')
else:
    print('Different distributions (reject H0)')<jupyter_output><empty_output><jupyter_text>2.3 Analysis of Variance (ANOVA)The same purpose as the T-Test but for more data samples.**Hypothesis**- $H_0$ = the mean across two or more groups are equal - **Fail to Reject $H_0$**: All sample distributions are equal - **Reject $H_0$**: One or more sample distributions are not equalURL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html<jupyter_code>from scipy.stats import f_oneway

seed(1)
alpha = 0.05

data1 = 5 * randn(100) + 50
data2 = 5 * randn(100) + 50
data3 = 5 * randn(100) + 52

# sns.distplot(data1, bins=10)
# sns.distplot(data2, bins=10)
# sns.distplot(data3, bins=10)

sns.histplot(data1, bins=10, color="blue")
sns.histplot(data2, bins=10, color="green")
sns.histplot(data3, bins=10, color="orange")

# compare samples
stat, p = f_oneway(data1, data2, data3)
print('Statistics=%.3f, p=%.3f' % (stat, p))

# interpret
if p > alpha:
    print('Same distributions (fail to reject H0)')
else:
    print('Different distributions (reject H0)')<jupyter_output><empty_output><jupyter_text>2.4 Mann-Whitney U testSimilar to T-test, but non-parametric for at least 20 observations in each data sample.**Hypothesis**- $H_0$ = no difference between the distributions of the data samples - **Fail to Reject $H_0$**: Sample distributions are equal - **Reject $H_0$**: Sample distributions are not equalURL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html<jupyter_code>from scipy.stats import mannwhitneyu

seed(1)
data1 = 10 * rand(100) + 50
data2 = 10 * rand(100) + 51

# sns.distplot(data1, bins=10)
# sns.distplot(data2, bins=10)

sns.histplot(data1, bins=10, color="blue")
sns.histplot(data2, bins=10, color="green")

# compare samples
stat, p = mannwhitneyu(data1, data2) 
print('Statistics=%.3f, p=%.3f' % (stat, p))

# interpret
alpha = 0.05
if p > alpha:
    print('Same distribution (fail to reject H0)') 
else:
    print('Different distribution (reject H0)')<jupyter_output><empty_output><jupyter_text>2.5 Kruskal-Wallis H TestThe Kruskal-Wallis test is a nonparametric version of the one-way analysis of variance test or ANOVA for short. This test can be used to determine whether more than two independent samples have a different distribution. It can be thought of as the generalization of the Mann-Whitney U test.**Hypothesis**- $H_0$ = all data samples were drawn from the same distribution - **Fail to Reject $H_0$**: All sample distributions are equal - **Reject $H_0$**: One or more sample distributions are not equal.URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kruskal.html<jupyter_code>from scipy.stats import kruskal

seed(1)
data1 = 5 * randn(100) + 50
data2 = 5 * randn(100) + 50
data3 = 5 * randn(100) + 52

# sns.distplot(data1, bins=10)
# sns.distplot(data2, bins=10)
# sns.distplot(data3, bins=10)

sns.histplot(data1, bins=10, color="blue")
sns.histplot(data2, bins=10, color="green")
sns.histplot(data3, bins=10, color="orange")

# compare samples
stat, p = kruskal(data1, data2, data3)
print('Statistics=%.3f, p=%.3f' % (stat, p))

# interpret
alpha = 0.05
if p > alpha:
    print('Same distributions (fail to reject H0)')
else:
    print('Different distributions (reject H0)')<jupyter_output><empty_output><jupyter_text>3. Statistical power analysis 3.1 Pearson's correlation between two samples $r_{xy} = \frac{ \sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y}) }{%        \sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2}\sqrt{\sum_{i=1}^{n}(y_i-\bar{y})^2}}$**Correlation**- value in the range ****- **Positive** correlation: both variables change in the same direction.- **Neutral** correlation: no relationship in the change of the variables. - **Negative** correlation: variables change in opposite directions.URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html<jupyter_code>from scipy.stats import pearsonr

seed(1)
data1 = 10 * randn(10000) + 50
data2 = data1 + (10 * randn(10000) + 50)

# plotting
# plt.plot(data1)
# plt.plot(data2)
sns.histplot(data1, bins=10, color="blue")
sns.histplot(data2, bins=10, color="green")

# calculate pearson's correlation
corr, _ = pearsonr(data1, data2) 
print('Pearsons correlation: %.3f' % corr)<jupyter_output><empty_output><jupyter_text>3.2 CovarianceCovariance provides the a measure of strength of correlation between two variable or more set of variables. The covariance matrix element $C_{ij}$ is the covariance of $x_i$ and $x_j$. The element $C_{ii}$ is the variance of $x_i$.- If $cov(x_i, x_j) = 0$ then variables are uncorrelated- If $cov(x_i, x_j) > 0$ then variables positively correlated- If $cov(x_i, x_j) < 0$ then variables negatively correlated$ np.cov(a,b) =    \begin{bmatrix}         cov(a,a) & cov(a,b) \\        cov(a,b) & cov(b,b) \\    \end{bmatrix}\quad$URL https://numpy.org/doc/stable/reference/generated/numpy.cov.html<jupyter_code>from numpy import cov

seed(1)
data1 = 20 * randn(1000) + 100
data2 = data1 + (10 * randn(1000) + 50)

# plotting
# sns.distplot(data1, bins=10)
# sns.distplot(data2, bins=10)

sns.histplot(data1, bins=10, color="blue")
sns.histplot(data2, bins=10, color="green")

# calculate covariance matrix
covariance = cov(data1, data2)[0, 1]
print(covariance)<jupyter_output><empty_output><jupyter_text>3.3 Effect size by Cohen's d $d = \frac{\mu_1 - \mu_2}{s}$using pooled standard deviation $s$ and variance $s^2$ $s = \sqrt{\frac{(n_1 - 1) s^2_{X_1} + (n_2 - 1) s^2_{X_2}}{n_1 + n_2  - 2}}$Cohen’s $d$ measures the difference between the mean from two Gaussian-distributed variables. Because the score is standardized, there is a table for the interpretation of the result, summarized as- Small effect: d=0.20- Medium effect: d=0.50 - **Large effect: d=0.80**<jupyter_code># calculate cohen's d for independent samples
def cohend(d1, d2):
    # calculate the size of samples
    n1, n2 = len(d1), len(d2) 
    
    # calculate the variance of the samples
    s1, s2 = var(d1, ddof=1), var(d2, ddof=1)

    # calculate the pooled standard deviation
    s = sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2)) 
    
    # calculate the means of the samples
    u1, u2 = mean(d1), mean(d2)
    
    # calculate the effect size
    d = (u1 - u2) / s 
    return d

# generate two samples
seed(1)
data1 = 5 * randn(100) + 52
data2 = 5 * randn(100) + 50

# data1 = stats.norm.rvs(loc=5, scale=10, size=500)
# data2 = stats.norm.rvs(loc=5, scale=10, size=500)

# plt.plot(data1)
# plt.plot(data2)
sns.distplot(data1, bins=10)
sns.distplot(data2, bins=10)

cd = cohend(data1, data2)
print('Cohens d value: %f' % cd)

# interpret
if 0.2 <= cd < 0.5:
    print('Small effect - Cohens d value: %f' % cd)
elif 0.5 <= cd < 0.8:
    print('Medium effect - Cohens d value: %f' % cd)
elif 0.8 <= cd:
    print('Large effect - Cohens d value: %f' % cd)<jupyter_output><empty_output><jupyter_text>3.4 Confidence interval if (n >= 30) then $\mu \pm z \frac{\sigma}{\sqrt{n}}$ if (n < 30)  then $\mu \pm t \frac{\sigma}{\sqrt{n}}$- Sample size $n$- Mean $\mu$- Standard deviation $\sigma$- Degrees of Freedom in a sample $df=n-1$- The z-value for 95% confidence is $Z=1.96$- Z-table https://www.statisticshowto.com/tables/z-table/    - t-table  https://www.statisticshowto.com/tables/t-distribution-table/URL https://www.statsmodels.org/devel/generated/statsmodels.stats.weightstats.DescrStatsW.tconfint_mean.html<jupyter_code>import statsmodels.stats.api as sms

print('data1: confidence interval', sms.DescrStatsW(data1).tconfint_mean())
print('data2: confidence interval', sms.DescrStatsW(data2).tconfint_mean())<jupyter_output><empty_output><jupyter_text>Combinated example of generated datasets**Populations**<jupyter_code>loc1 = 80
loc2 = 78
s = 5
data1 = stats.norm(loc=loc1, scale=s)
data2 = stats.norm(loc=loc2, scale=s)

start = 60
stop = 100
num = 1000
x = np.linspace(start, stop, num)
plt.plot(x, data1.pdf(x), 'b')
plt.plot(x, data2.pdf(x), 'r')<jupyter_output><empty_output><jupyter_text>**Samples**<jupyter_code>sample_size = 100
sample1 = data1.rvs(sample_size)
sample2 = data2.rvs(sample_size)
print('sample1 mean=', sample1.mean())
print('sample2 mean=', sample2.mean())
sns.distplot(sample1, bins=5)
sns.distplot(sample2, bins=5)<jupyter_output><empty_output><jupyter_text>**Test of normality: if samples have nornal distribution**<jupyter_code>print('Shapiro sample1:', stats.shapiro(sample1))
print('Shapiro sample2:', stats.shapiro(sample1))<jupyter_output><empty_output><jupyter_text>**Test of variances: if samples are from populations with equal variances**<jupyter_code>stats.levene(sample1, sample2)<jupyter_output><empty_output><jupyter_text>**Statistical power**<jupyter_code>cd = cohend(sample1, sample2)
print('Cohens d value: %f' % cd)<jupyter_output><empty_output><jupyter_text>**Estimate optimal sample size based on the Cohen's d value**- solve for any one parameter of the power of a two sample t-test- for t-test the keywords are: effect_size, nobs1, alpha, power, ratio- exactly one needs to be None, all others need numeric valuesURL https://www.statsmodels.org/stable/generated/statsmodels.stats.power.tt_ind_solve_power.html<jupyter_code>sm_stats.power.tt_ind_solve_power(cd, None, 0.05, 0.8, 1)<jupyter_output><empty_output><jupyter_text>**Two-sided confidence interval for weighted mean of data**<jupyter_code>print('sample1: confidence interval', sms.DescrStatsW(sample1).tconfint_mean())
print('sample2: confidence interval', sms.DescrStatsW(sample2).tconfint_mean())<jupyter_output><empty_output><jupyter_text>**Experiment toward reality ;)**<jupyter_code>df = pd.DataFrame({'observation': np.repeat([True, False], sample_size), 
                   'score': np.concatenate((sample1, sample2))})
sns.boxplot(x='observation', y='score', data=df)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14359279065161418, 'fraction_numerical': 0.04502541757443718, 'mean_word_length': 5.198036006546645, 'pattern_counts': {'":': 0, '<': 109, '<?xml version=': 0, '>': 109, 'https://': 14, 'lorem ipsum': 0, 'www.': 4, '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': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '154112', 'license': 'Apache-2.0', 'n_tokens_mistral': 6046, 'n_tokens_neox': 5562, 'n_words': 1444, 'path': 'cvicenia/tyzden-05/IAU_051_EDA_stats.ipynb', 'repo_name': 'FIIT-IAU/2021-2022', 'size': '15147'}