Распределение Пуассона в Python: методы расчета и визуализация

P(X = k) = (λ^k * e^(-λ)) / k!

import numpy as np
from scipy import stats

# Параметр распределения (среднее число событий)
lambda_param = 3.5

# Вычисление вероятности ровно 4 событий
k = 4
probability = stats.poisson.pmf(k, lambda_param)
print(f"Вероятность ровно {k} событий при λ={lambda_param}: {probability:.6f}")

# Вычисление вероятности не более 2 событий (0, 1 или 2)
probability_cumulative = stats.poisson.cdf(2, lambda_param)
print(f"Вероятность не более 2 событий при λ={lambda_param}: {probability_cumulative:.6f}")

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# Генерация 1000 случайных наблюдений из распределения Пуассона с λ=5
np.random.seed(42) # Для воспроизводимости результатов
sample = stats.poisson.rvs(mu=5, size=1000)

# Оценка параметра λ из выборки
lambda_estimate = np.mean(sample)
print(f"Истинное значение λ: 5")
print(f"Оценка λ из выборки: {lambda_estimate:.4f}")

# Подгонка распределения Пуассона к данным
lambda_fit = stats.poisson.fit(sample)
print(f"Оценка λ методом максимального правдоподобия: {lambda_fit[0]:.4f}")

import numpy as np
from scipy import stats

lambda_param = 4.2
# Вероятность от 3 до 7 событий включительно
probability_range = stats.poisson.pmf(np.arange(3, 8), lambda_param).sum()
print(f"Вероятность от 3 до 7 событий при λ={lambda_param}: {probability_range:.6f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Диапазон возможных значений k
k_values = np.arange(0, 20)

# Различные значения параметра λ
lambda_values = [1\.5, 3, 6, 10]

plt.figure(figsize=(12, 8))

for lambda_param in lambda_values:
# Вычисление вероятностей для каждого k
pmf_values = stats.poisson.pmf(k_values, lambda_param)

# Построение графика
plt.plot(k_values, pmf_values, 'o-', label=f'λ = {lambda_param}')

plt.title('Распределение Пуассона для различных значений λ', fontsize=15)
plt.xlabel('Число событий (k)', fontsize=12)
plt.ylabel('Вероятность P(X = k)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Параметр распределения
lambda_param = 5

# Генерация выборки из распределения Пуассона
np.random.seed(42)
sample = stats.poisson.rvs(lambda_param, size=1000)

# Теоретические вероятности
k_values = np.arange(0, 15)
pmf_theoretical = stats.poisson.pmf(k_values, lambda_param)

plt.figure(figsize=(12, 7))

# Гистограмма эмпирических данных
plt.hist(sample, bins=np.arange(0, 15+1)-0.5, density=True, alpha=0.5, label='Эмпирическое распределение')

# Теоретическое распределение
plt.plot(k_values, pmf_theoretical, 'ro-', label='Теоретическое распределение')

plt.title(f'Сравнение теоретического и эмпирического распределений Пуассона (λ={lambda_param})', fontsize=14)
plt.xlabel('Число событий (k)', fontsize=12)
plt.ylabel('Вероятность / Частота', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=12)
plt.xticks(k_values)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

lambda_param = 4

# Диапазон значений k
k_values = np.arange(0, 15)

# Создание сетки из двух графиков
fig, axes = plt.subplots(2, 1, figsize=(10, 12))

# PMF – функция вероятности массы
pmf_values = stats.poisson.pmf(k_values, lambda_param)
axes[0].bar(k_values, pmf_values, alpha=0.7, width=0.8)
axes[0].set_title(f'PMF распределения Пуассона (λ={lambda_param})', fontsize=14)
axes[0].set_xlabel('Число событий (k)', fontsize=12)
axes[0].set_ylabel('Вероятность P(X = k)', fontsize=12)
axes[0].grid(True, alpha=0.3)

# CDF – кумулятивная функция распределения
cdf_values = stats.poisson.cdf(k_values, lambda_param)
axes[1].step(k_values, cdf_values, 'r-', where='post', linewidth=2)
axes[1].scatter(k_values, cdf_values, c='r', s=50)
axes[1].set_title(f'CDF распределения Пуассона (λ={lambda_param})', fontsize=14)
axes[1].set_xlabel('Число событий (k)', fontsize=12)
axes[1].set_ylabel('Вероятность P(X ≤ k)', fontsize=12)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Создание сетки графиков 2x2
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

# Различные значения λ для визуализации
lambda_values = [1, 5, 10, 20]

for i, lambda_param in enumerate(lambda_values):
k_range = np.arange(0, max(lambda_param * 2, 10))
pmf_values = stats.poisson.pmf(k_range, lambda_param)

axes[i].bar(k_range, pmf_values, alpha=0.7)
axes[i].set_title(f'λ = {lambda_param}', fontsize=14)
axes[i].set_xlabel('Число событий (k)', fontsize=12)
axes[i].set_ylabel('Вероятность P(X = k)', fontsize=12)
axes[i].grid(True, alpha=0.3)

# Отображение среднего (λ) на графике
axes[i].axvline(lambda_param, color='red', linestyle='--', 
label=f'Среднее (λ={lambda_param})')
axes[i].legend(fontsize=10)

plt.suptitle('Форма распределения Пуассона при различных значениях λ', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Среднее количество звонков в час
lambda_calls = 15

# Диапазон для анализа
k_range = np.arange(0, 30)

# Расчет вероятностей
pmf_values = stats.poisson.pmf(k_range, lambda_calls)

# Вероятность более 20 звонков в час
prob_more_than_20 = stats.poisson.sf(20, lambda_calls)
print(f"Вероятность более 20 звонков в час: {prob_more_than_20:.4f} ({prob_more_than_20*100:.2f}%)")

# Минимальное количество операторов для обработки звонков с 95% уверенностью
min_operators = stats.poisson.ppf(0.95, lambda_calls)
print(f"Минимальное число операторов (95% уверенность): {np.ceil(min_operators)}")

# Визуализация
plt.figure(figsize=(10, 6))
plt.bar(k_range, pmf_values, alpha=0.7)
plt.axvline(min_operators, color='r', linestyle='--', 
label=f'95% уверенность: {np.ceil(min_operators)} операторов')
plt.title('Распределение числа звонков в час', fontsize=14)
plt.xlabel('Количество звонков', fontsize=12)
plt.ylabel('Вероятность', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Среднее число дефектов на 100 м
lambda_defects = 3.5

# Симуляция контроля 50 рулонов ткани по 100 м
np.random.seed(42)
defects_per_roll = stats.poisson.rvs(lambda_defects, size=50)

# Расчет контрольных пределов (на основе 3-сигма)
ucl = lambda_defects + 3 * np.sqrt(lambda_defects) # Верхний контрольный предел
lcl = max(0, lambda_defects – 3 * np.sqrt(lambda_defects)) # Нижний контрольный предел

# Вероятность выхода за UCL
prob_above_ucl = stats.poisson.sf(int(ucl), lambda_defects)
print(f"Верхний контрольный предел: {ucl:.2f} дефектов")
print(f"Вероятность превышения UCL: {prob_above_ucl:.6f} ({prob_above_ucl*100:.4f}%)")

# Построение контрольной карты
plt.figure(figsize=(12, 6))
plt.plot(range(1, 51), defects_per_roll, 'bo-', alpha=0.7)
plt.axhline(lambda_defects, color='g', linestyle='-', label=f'Среднее (λ={lambda_defects})')
plt.axhline(ucl, color='r', linestyle='--', label=f'UCL={ucl:.2f}')
plt.axhline(lcl, color='r', linestyle='--', label=f'LCL={lcl:.2f}')

# Подсветка точек вне контрольных пределов
out_of_control = np.where((defects_per_roll > ucl) | (defects_per_roll < lcl))[0]
if len(out_of_control) > 0:
plt.plot(out_of_control + 1, defects_per_roll[out_of_control], 'ro', markersize=10)

plt.title('Контрольная карта дефектов по рулонам ткани', fontsize=14)
plt.xlabel('Номер рулона', fontsize=12)
plt.ylabel('Количество дефектов', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Среднее количество запросов в минуту
lambda_requests = 120

# Оценка вероятности пиковой нагрузки (>150 запросов/мин)
prob_peak = stats.poisson.sf(150, lambda_requests)
print(f"Вероятность более 150 запросов в минуту: {prob_peak:.4f} ({prob_peak*100:.2f}%)")

# Симуляция трафика на протяжении часа (60 минут)
np.random.seed(42)
traffic_per_minute = stats.poisson.rvs(lambda_requests, size=60)

# Визуализация
plt.figure(figsize=(12, 7))
plt.bar(range(1, 61), traffic_per_minute, alpha=0.7)
plt.axhline(lambda_requests, color='g', linestyle='-', 
label=f'Среднее (λ={lambda_requests})')
plt.axhline(150, color='r', linestyle='--', 
label=f'Порог пиковой нагрузки (150 запросов)')

plt.title('Моделирование трафика веб-сайта', fontsize=14)
plt.xlabel('Минута', fontsize=12)
plt.ylabel('Количество запросов', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats
import math
import time

def poisson_pmf_scipy(k, lambda_param):
"""Расчет через scipy.stats."""
return stats.poisson.pmf(k, lambda_param)

def poisson_pmf_numpy(k, lambda_param):
"""Расчет через numpy."""
return np.exp(-lambda_param) * np.power(lambda_param, k) / np.math.factorial(k)

def poisson_pmf_math(k, lambda_param):
"""Расчет через стандартную библиотеку math."""
return math.exp(-lambda_param) * (lambda_param ** k) / math.factorial(k)

def poisson_pmf_log(k, lambda_param):
"""Расчет с использованием логарифмов для предотвращения переполнения."""
log_p = k * np.log(lambda_param) – lambda_param – np.log(math.factorial(k))
return np.exp(log_p)

def poisson_pmf_vectorized(k_array, lambda_param):
"""Векторизованный расчет для массива k."""
log_p = k_array * np.log(lambda_param) – lambda_param – np.array([np.log(math.factorial(k)) for k in k_array])
return np.exp(log_p)

# Тестирование производительности для одного значения k
k = 10
lambda_param = 5
repetitions = 100000

methods = {
'scipy': lambda: poisson_pmf_scipy(k, lambda_param),
'numpy': lambda: poisson_pmf_numpy(k, lambda_param),
'math': lambda: poisson_pmf_math(k, lambda_param),
'log': lambda: poisson_pmf_log(k, lambda_param)
}

results = {}

for name, method in methods.items():
start_time = time.time()
for _ in range(repetitions):
method()
results[name] = time.time() – start_time

print("Время выполнения для одиночного расчета (секунды):")
for name, elapsed in results.items():
print(f"{name}: {elapsed:.6f}")

# Тестирование для массива значений k
k_array = np.arange(0, 100)
repetitions_array = 1000

methods_array = {
'scipy_loop': lambda: [stats.poisson.pmf(k, lambda_param) for k in k_array],
'scipy_vectorized': lambda: stats.poisson.pmf(k_array, lambda_param),
'custom_vectorized': lambda: poisson_pmf_vectorized(k_array, lambda_param)
}

results_array = {}

for name, method in methods_array.items():
start_time = time.time()
for _ in range(repetitions_array):
method()
results_array[name] = time.time() – start_time

print("\nВремя выполнения для массива расчетов (секунды):")
for name, elapsed in results_array.items():
print(f"{name}: {elapsed:.6f}")

import numpy as np
from scipy import stats
import time
import matplotlib.pyplot as plt

# Большой диапазон значений λ
lambda_values = np.linspace(0.1, 50, 1000)
k = 10 # Интересуемся вероятностью ровно 10 событий

# Неоптимальный подход (создает большой промежуточный массив)
def calculate_unoptimized():
return stats.poisson.pmf(k, lambda_values)

# Оптимизированный подход (использует генератор)
def calculate_optimized():
return np.array([stats.poisson.pmf(k, lam) for lam in lambda_values])

# Сравнение производительности
print("Измерение производительности для %d значений λ:" % len(lambda_values))

start_time = time.time()
result1 = calculate_unoptimized()
unopt_time = time.time() – start_time
print(f"Неоптимизированный метод: {unopt_time:.6f} сек")

start_time = time.time()
result2 = calculate_optimized()
opt_time = time.time() – start_time
print(f"Оптимизированный метод: {opt_time:.6f} сек")

print(f"Разница составляет {unopt_time/opt_time:.2f}x")

# Визуализация результатов
plt.figure(figsize=(12, 6))
plt.plot(lambda_values, result1, label='PMF для k=10')
plt.axvline(k, color='r', linestyle='--', label='λ = k (максимум PMF)')
plt.title('Вероятность ровно 10 событий в зависимости от λ', fontsize=14)
plt.xlabel('Параметр λ', fontsize=12)
plt.ylabel('Вероятность P(X = 10)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# Большое значение λ
lambda_param = 1000

# Диапазон значений k для визуализации
k_range = np.arange(lambda_param – 4 * np.sqrt(lambda_param), 
lambda_param + 4 * np.sqrt(lambda_param), dtype=int)

# Точное вычисление через распределение Пуассона
poisson_pmf = stats.poisson.pmf(k_range, lambda_param)

# Аппроксимация нормальным распределением
normal_approx = stats.norm.pdf(k_range, loc=lambda_param, scale=np.sqrt(lambda_param))

plt.figure(figsize=(12, 6))
plt.plot(k_range, poisson_pmf, 'ro-', alpha=0.7, label='Точное распределение Пуассона')
plt.plot(k_range, normal_approx, 'b-', label='Нормальная аппроксимация')
plt.title(f'Аппроксимация распределения Пуассона (λ={lambda_param}) нормальным распределением', fontsize=14)
plt.xlabel('Число событий (k)', fontsize=12)
plt.ylabel('Вероятность P(X = k)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from concurrent.futures import ProcessPoolExecutor
import time

# Функция для вычисления вероятностей для одного λ
def calculate_stats_for_lambda(lambda_param):
k_range = np.arange(0, int(lambda_param * 2) + 1)
pmf = stats.poisson.pmf(k_range, lambda_param)
mean = np.sum(k_range * pmf)
variance = np.sum((k_range – mean)**2 * pmf)
return lambda_param, mean, variance, np.max(pmf)

# Последовательный расчет
def sequential_calculation(lambda_values):
results = []
for lam in lambda_values:
results.append(calculate_stats_for_lambda(lam))
return results

# Параллельный расчет
def parallel_calculation(lambda_values, max_workers=None):
with ProcessPoolExecutor(max_workers=max_workers) as executor:
results = list(executor.map(calculate_stats_for_lambda, lambda_values))
return results

# Тестирование производительности
lambda_values = np.linspace(1, 100, 100)

print("Последовательный расчет:")
start_time = time.time()
seq_results = sequential_calculation(lambda_values)
seq_time = time.time() – start_time
print(f"Время выполнения: {seq_time:.6f} сек")

print("\nПараллельный расчет:")
start_time = time.time()
par_results = parallel_calculation(lambda_values)
par_time = time.time() – start_time
print(f"Время выполнения: {par_time:.6f} сек")
print(f"Ускорение: {seq_time/par_time:.2f}x")