Статистика с Python в розничной торговле : нормальное приближение и биномиальное распределение

Задачи этого модуля :

• Эмпирическое правило нормального распределения
• Знакомство с биномиальным распределением и его нормальной аппроксимацией
• Стандартизация данных и нормальное приближение
• Центральная предельная теорема

Вернемся к гистограмме длины чека,  это типичный пример колоколообразной гистограммы. В таком случае говорят, что данные следуют нормальному закону распределения. Однако далеко не все гистограммы выглядят подобным образом, мы это видели на примере распределения выручек магазинов.

Загружаем данные

df=pd.read_excel("shops_lch_07.xlsx",index_col="shop")

df

lch
shop
Sh1	3.617
Sh2	3.513
...	...
Sh127	3.045
Sh128	3.196

128 rows × 1 columns

Выводим гистограмму с обозначением диапазонов стандартного отклонения

plt.figure(figsize=(10,6))
mean = df.lch.mean()
median =df.lch.median()
sd_1=df.lch.std()
sd_2=2*df.lch.std()
sd_3=3*df.lch.std()
sns.histplot(df,stat="percent",kde=True)
plt.axvline(x=mean, c='g', label='Mean')
plt.axvline(x=median, c='black', label='Median')
plt.axvline(x=mean-sd_1, c='red',label='mean-1*sd',ls='-')
plt.axvline(x=mean+sd_1, c='red',label='mean+1*sd',ls='-')
plt.axvline(x=mean-sd_2, c='red',label='mean-2*sd',ls='--')
plt.axvline(x=mean+sd_2, c='red',label='mean+2*sd',ls='--')
plt.axvline(x=mean-sd_3, c='red',label='mean-3*sd',ls='-.')
plt.axvline(x=mean+sd_3, c='red',label='mean+3*sd',ls='-.')
plt.legend()

Эмпирическое правило

Поскольку все нормальные распределения имеют одинаковую общую форму, мы можем сформулировать некоторые суждения о том, как распределены данные при любом нормальном распределении. Эмпирическое правило гласит, что для любого нормального распределения:

• около 68% данных находятся в интервале ± одно стандартное отклонение от среднего;

• около 95% данных находятся в интервале ± два стандартных отклонения от среднего;

• около 99% данных находятся в интервале ± три стандартных отклонения от среднего.

Попробуем проверить это на наших данных

pct_1sd=df[(df.lch>(mean-1*sd_1))&(df.lch<(mean+1*sd_1))].count()[0]/df.lch.count()

pct_2sd=df[(df.lch>(mean-2*sd_1))&(df.lch<(mean+2*sd_1))].count()[0]/df.lch.count()

pct_3sd=df[(df.lch>(mean-3*sd_1))&(df.lch<(mean+3*sd_1))].count()[0]/df.lch.count()

print(f'percentage of data is within ± 1 standard deviation of the mean = {pct_1sd:.2f}')

print(f'percentage of data is within ± 2 standard deviation of the mean = {pct_2sd:.2f}')

print(f'percentage of data is within ± 3 standard deviation of the mean = {pct_3sd:.2f}')

percentage of data is within ± 1 standard deviation of the mean = 0.84
percentage of data is within ± 2 standard deviation of the mean = 0.98
percentage of data is within ± 3 standard deviation of the mean = 0.98

Как видно, на наших данных эмпирическое правило не работает,

то есть наши данные не соответствуют нормальному распределению и причина этому - 

наличие выбросов. Устраним их. Для этого воспользуемся правилом трех стандартных

отклонений

Три стандартных отклонения от среднего значения - это обычное ограничение на практике для выявления выбросов при нормальном распределении. Если распределение скошено, то для определения выбросов лучше использовать правило межквартильного размаха.

Убираем выбросы и выводим новую гистограмму

norm = (df.lch < (mean+sd_3)) & (df.lch> (mean-sd_3))

df_norm=df[norm]

mean_n = df_norm.lch.mean()

median_n =df_norm.lch.median()

sdn_1=df_norm.lch.std()

sdn_2=2*df_norm.lch.std()

sdn_3=3*df_norm.lch.std()

plt.figure(figsize=(10,6))

sns.histplot(df_norm,stat="percent",kde=True)

plt.axvline(x=mean_n, c='g', label='Mean')

plt.axvline(x=median_n, c='black', label='Median')

plt.axvline(x=mean_n-sdn_1, c='red',label='mean-1*sd',ls='-')

plt.axvline(x=mean_n+sdn_1, c='red',label='mean+1*sd',ls='-')

plt.axvline(x=mean_n-sdn_2, c='red',label='mean-2*sd',ls='--')

plt.axvline(x=mean_n+sdn_2, c='red',label='mean+2*sd',ls='--')

plt.axvline(x=mean_n-sdn_3, c='red',label='mean-3*sd',ls='-.')

plt.axvline(x=mean_n+sdn_3, c='red',label='mean+3*sd',ls='-.')

plt.legend()

На этих данных эмпирическое правило работает

pct_1sd=df_norm[(df_norm.lch>(mean_n-1*sdn_1))&(df_norm.lch<(mean_n+1*sdn_1))].count()[0]/df_norm.lch.count()

pct_2sd=df_norm[(df_norm.lch>(mean_n-2*sdn_1))&(df_norm.lch<(mean_n+2*sdn_1))].count()[0]/df_norm.lch.count()

pct_3sd=df_norm[(df_norm.lch>(mean_n-3*sdn_1))&(df_norm.lch<(mean_n+3*sdn_1))].count()[0]/df_norm.lch.count()

print(f'percentage of data is within ± 1 standard deviation of the mean = {pct_1sd:.2f}')

print(f'percentage of data is within ± 2 standard deviation of the mean = {pct_2sd:.2f}')

print(f'percentage of data is within ± 3 standard deviation of the mean = {pct_3sd:.2f}')

percentage of data is within ± 1 standard deviation of the mean = 0.70
percentage of data is within ± 2 standard deviation of the mean = 0.94
percentage of data is within ± 3 standard deviation of the mean = 1.00

Стандартизация данных и стандартная нормальная кривая
Стандартное нормальное распределение — это такое распределение, в котором
единицы на оси х выражены в стандартных отклонениях от среднего. Для того чтобы
сравнить данные со стандартным нормальным распределением, нужно вычесть
среднее и затем разделить на стандартное отклонение; эта процедура также называется
нормализацией или стандартизацией. Преобразованное значение называется

z-оценкой, или стандартной оценкой, а нормальное распределение иногда называют

z-распределением.
С помощью стандартизации отбрасывание (или отчистка) выбросов наших данных

 могла бы  выглядеть следующим образом

Для стандартизации используем функцию zscore из пакета scipy

from scipy import stats

z = np.abs(stats.zscore(df))

clean_df = df

clean_df = clean_df[z.lch < 3]

Сравниваем два полученных DataFrame

clean_df.lch.compare(df_norm.lch,keep_shape=True,keep_equal=True)

self	other
shop
Sh1	3.617	3.617
Sh2	3.513	3.513
...	...	...
Sh127	3.045	3.045
Sh128	3.196	3.196

126 rows × 2 columns

Нормальное приближение

Для начала введем понятие кумулятивная функция распределения

(cumulative distribution function, CDF), которая дает вероятность, что случайная величина

меньше или равна некоторому значению. Для нашего случая с длиной чека вероятность равносильна доле чеков с определенной длиной. 

Нормальная аппроксимация означает использование площади под нормальной кривой для вычисления процентов. Давайте посмотрим на пример. Какой процент магазинов имеют длину чека от 3 до 3.5 ?Помните, что мы знаем среднее значение и стандартное отклонение для этих данных, и это действительно все, что нам нужно знать.

Создаем экземпляр объекта случайной величины с нормальным распределением и с параметрами наших данных, нанесем ее график на нашу гистограмму и обозначим диапазон, площадь которую хотим определить

lch_rv = stats.norm(loc = mean_n,scale = sdn_1)

lch_rv_pdf = lch_rv.pdf(values)

plt.figure(figsize=(10,6))

values = np.linspace(2, 4, num=100)

sns.histplot(df_norm,stat="density")

plt.plot(values, lch_rv_pdf)

plt.axvline(x=3, c='red',label='lch=3',ls='-')

plt.axvline(x=3.5, c='red',label='lch=3.5',ls='-')

plt.grid()

plt.legend()

А теперь с помощью CDF функции рассчитаем площадь в процентах под нормальной кривой и сравним с истинным процентом наших данных в этом диапазоне

upper_lim=3.5

lower_lim=3

prob_lch_in_3_3_5 = lch_rv.cdf(upper_lim) - lch_rv.cdf(lower_lim)

pct_lch_in_3_3_5 = df_norm[(df_norm.lch>lower_lim)&(df_norm.lch<upper_lim)].count()[0]/df_norm.lch.count()

print(f'prob(3 <= lch <= 3.5) = {prob_lch_in_3_3_5:0.4f}')

print(f'pct(3 <= lch <= 3.5) = {pct_lch_in_3_3_5:0.4f}')

prob(3 <= lch <= 3.5) = 0.4665
pct(3 <= lch <= 3.5) = 0.4444

Как видим, нормальное приближение получилось довольно точное.

Вычисление перцентилей с помощью нормального приближения

Определим новое понятие - квантильная функция, также известная как обратная CDF
или процентная функция (PPF), связана с распределением случайных величин.

В то время как CDF сообщает нам вероятность того, что случайная величина X

должна иметь значение, равное или меньшее определенное значение x,

квантильная функция сообщает нам значение случайной величины такая, что

вероятность того, что переменная меньше или равно этому значению, равно

заданной вероятности.

Теперь посмотрим на другой пример использования нормального приближения.

Напомним, что 30-й перцентиль — это точка, ниже которой находится 30% данных. Если мы нарисуем нормальную кривую, тогда 30-й перцентиль будет значением, для 30%

данных попадают ниже. Найдем с помощью нормального приближения значение для

30-го перцентиля

print(f'30th percentile value function ppf = {lch_rv.ppf(0.30):0.4f}')
print(f'30th percentile value function quantile= \ 

{df_norm.lch.quantile(0.30):0.4f}')

30th percentile value function ppf = 2.8597
30th percentile value function quantile= 2.8977

Биномиальное распределение и нормальное приближение к биномиальному

Биномиальное распределение-это теоретическое распределение, которое применяется к случайным переменным и удовлетворяет следующим трем характеристикам:

• Условие 1: для индивидуального наблюдения существует только два возможных            результата, обычно обозначаемых как успех и неудача. Если вероятность успеха
равна p, то вероятность неудачи должна быть равна 1 – p.

• Условие 2: эксперимент выполняется фиксированное число раз, обычно
обозначаемое n.

• Условие 3: все эксперименты независимы, что означает, что знание результата             эксперимента не меняет вероятности следующего.

Поэтому вероятность успеха (и неудачи) остается прежней.
Если эти условия выполняются, то мы говорим, что случайная величина следует             биномиальному распределению  или что случайная величина является биномиальной   случайной величиной.

Пример работы с биномиальным распределением :

Компания электронной коммерции продает женскую одежду. Замечено, что
10% их клиентов возвращают товары, купленные ими по многим причинам (таким как несоответствие размера, цвета и материала). В определенный день 20 клиентов                приобрели товары у этой компании . 

Вычислить:

1. Вероятность того, что ровно 5 клиентов вернут товар.
2. Вероятность того, что максимум 5 клиентов вернут товар.
3. Вероятность того, что более 5 клиентов будут возвращать товары, приобретенные по ним.
4. Среднее число клиентов, которые, скорее всего, вернут товары, а также дисперсия и
стандартное отклонение числа возвратов.

Мы решаем каждый из них следующим образом:

1. Вероятность того, что ровно 5 клиентов вернут товар.

Функция stats.binom.pmf() вычисляет PMF для биномиального распределения и 
принимает три параметра:
а) ожидаемое число успешных испытаний (5)
b) общее число испытаний (20)
c) вероятность успеха (0,1)

from scipy import stats
stats.binom.pmf(5, 20, 0.1)
0.03192136111995428

Соответствующая вероятность равна 0,03192, то есть вероятность того, 
что ровно 5 клиентов вернут товар, составляет примерно 3%.
Чтобы визуализировать, как изменяется PMF с увеличением числа успешных
испытаний, мы создадим список всех возможных чисел успехов (от 0 до 20) и 
соответствующих значений PMF и нарисуем гистограмму, как показано на рис. 3.1.

# range(0,21) возвращает все значения от 0 до 20 (исключая 21)
pmf_df = pd.DataFrame({ 'success': range(0,21),'pmf': list(stats.binom.pmf(range(0,21),20, 0.1))})
# Выводим cтолбчатую диаграмму с количеством успехов и их вероятностью
plt.figure(figsize=(10,6))
sns.barplot(x = pmf_df.success, y = pmf_df.pmf)
plt.ylabel('pmf')
plt.xlabel('Number of items returned')



2. Вероятность того, что максимум 5 клиентов вернут товар.

Функция stats.binom.cdf() вычисляет CDF для биномиального распределения. В этом      случае функция кумулятивного распределения возвращает вероятность того,
что максимум 5 клиентов вернут товары.

stats.binom.cdf(5, 20, 0.1)
0.988746865835491

3. Вероятность того, что более 5 клиентов будут возвращать товары

Общая вероятность любого количества клиентов, возвращающих товары (включая 0),    всегда равна 1,0. таким образом, вероятность того, что более 5 клиентов вернут
товары, может быть вычислена путем вычитания вероятности того, что максимум 
5 клиентов  вернут товары из 1,0. Другими словами, вероятность того, что более
5 клиентов вернут  товары, может быть получена путем вычисления CDF из 5, а затем
вычитания его из 1.0.

1 - stats.binom.cdf(5, 20, 0.1)
0.011253134164509015

4. среднее число клиентов, которые, вероятно, вернут товары, а также дисперсия и
стандартное отклонение числа возвратов.

а) среднее значение биномиального распределения задается  n * p
(b) дисперсия биномиального распределения задается  n * p * (1 − p)

mean, var = stats.binom.stats(20, 0.1)
print('Average: ', mean , 'Variance:', var)
Average:  2.0 Variance: 1.8

Нормальное распределение также можно рассматривать как непрерывный предел
биномиального распределения n→∞. Подтверждение этому можно увидеть это на графике:

from scipy.stats import binom
cols = colors.cnames
n_values = [1, 5,10, 30, 100]
subplots = [111+100*x for x in range(0,len(n_values))]
ctr = 0
fig, ax = plt.subplots(len(subplots), figsize=(6,12))
k = np.arange(0, 200)
p=0.5
for n in n_values:
    k=np.arange(0,n+1)
    rv = binom(n, p)
    ax[ctr].plot(k, rv.pmf(k), lw=2)
    ax[ctr].set_title("n=" + str(n))
    ctr += 1
    fig.subplots_adjust(hspace=0.5)
    plt.suptitle("Binomial distribution PMF (p=0.5) for various values of n", fontsize=14)



Центральная предельная теорема (central limit theorem)
В соответствии с центральной предельной теоремой, при увеличении размера выборки
средние по выборке начинают подчиняться нормальному закону распределения и
концентрироваться вокруг среднего по совокупности. Это утверждение оказывается
верно независимо от распределения совокупности, из которой была получена
выборка. Вернемся к распределению выручки 128 магазинов, как мы видели ранее,
распределение скошено и далеко не похоже на нормальное. Помимо  самих данных
возьмем из них  две выборки по 128 значений  средних по 5 значений
и 128 средних по 20 значений. Затем построим три гистограммы.

sample_data = pd.DataFrame({
'rev': df.rev,
'type': 'fact',
})
sample_mean_05 = pd.DataFrame({
'rev': [df.rev.sample(5).mean() for _ in range(128)],
'type': 'average sample of 5',
})
sample_mean_20 = pd.DataFrame({
'rev': [df.rev.sample(20).mean() for _ in range(128)],
'type': 'average sample of 20',
})
results = pd.concat([sample_data, sample_mean_05, sample_mean_20])
plt.figure(figsize=(10,6))
g = sns.FacetGrid(results, col='type', col_wrap=1, height=4, aspect=2)
#g.map(plt.hist, 'rev',  bins=20)
g.map(sns.histplot, 'rev')
g.map(plt.axvline,x=df.rev.mean(), c='red',ls='--')
g.set_axis_labels('revenue', 'count')
g.set_titles('{col_name}')


Получилась хорошая иллюстрация центральной предельной теоремы.