ЛР / ЛР3 / SM_LR3

from matplotlib import pyplot as plt
import numpy as np


def random_erlang_val(k, lmbda, N=10000):
    '''Функция создающая случайное число распределенное по закону Эрланговского распределения'''
    f = lambda : -1/lmbda * sum(np.log(np.random.uniform(0,1)) for i in range(k))
    return [f() for n in range(N)]


def slicearray(array, intervals_cnt = 10):
    '''Функция для разбивки списка на интервалы и подсчета кол. попаданий'''
    interval=[] # Список интервалов
    entries = [] # Список кол-ва вхождений
    x_min, x_max = min(array), max(array)
    delta = (x_max - x_min) / intervals_cnt
    
    # Заполняем список интервалов
    for x in range (intervals_cnt+1):
        interval.append(x_min + delta * x)

    # Заполняем список кол-ва вхождений
    for i in range(intervals_cnt):
        entries.append(sum(True for x in array if interval[i] < x <= interval[i+1]))
    return interval, entries, x_min, x_max, delta


def ploting_erlang_distribution(k, lmbda, data = [], N=10000, iner_cnt = 50, add_hist=False):
    '''Построение графика плотности Эрланговского распределения'''
    plt.figure(figsize=(9,7))

    # Проверка переданны ли дынные для построения гистограммы
    if not data: data = random_erlang_val(k, lmbda, N)
    else: N = len(data)

    # Определение интервалов и их длительности, min/max значений для построения гистограммы
    intervals, entries, x_min, x_max, delta = slicearray(data, iner_cnt)

    # Проверка добавлять гистограмму или нет
    if add_hist:
        # Подсчет весов столбцов гистограммы
        weights = [cnt / (N * delta) for cnt in entries]
        plt.hist(intervals[:-1], intervals, weights = weights, label='Рассчитанное значение')

    # Функция плотности вероятности эрланговского распределения для заданных k и lambda
    f = lambda x : ((lmbda**k) * (x**(k-1)) * np.exp(-lmbda*x)) / np.math.factorial(k-1)

    # Генерируем список значений x для построения графика
    x = np.linspace(0 if x_min<=1 else x_min-1, x_max+1, 100)
    plt.plot(x, f(x), label='График плотности вероятности\nЭрланговского распределения\nдля k={}, lambda={}'.format(k, lmbda), color='red')
    plt.title("График плотности вероятности\nЭрланговского распределения", fontsize=16)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.legend()
    plt.show()


def ploting_lambda_nu(k, lmbda, lmbda_list, nu_list, N_array):
    '''Ф-ция построения графиков зависимости оценок интенсивности и коэффициента вариации от объема выборки заданного потока'''
    # Расчет теоретических значений
    lambda_t = lmbda / k
    nu_t = 1 / np.sqrt(k)

    # Создаем объекты графиков
    fig, axis = plt.subplots(1, 2, figsize=(14, 7))
    fig.suptitle('Статистика для Эрланговского распределения с k={} и lambda={}'.format(k, lmbda), fontsize=16)

    # Строим график зависимости оценок интенсивности от объема выборки
    axis[0].axhline(lambda_t, color='red', label='Теоретическое значение')
    axis[0].plot(N_array, lmbda_list, label='Рассчитанное значение')
    axis[0].set_title('График зависимости\nоценок интенсивности от объема выборки')
    axis[0].set_xlabel('Размер выборки')
    axis[0].set_ylabel('Оценка интенсивности')
    axis[0].legend()

    # Строим график зависимости оценок коэффициента вариации от объема выборки
    axis[1].axhline(nu_t, color='red', label='Теоретическое значение')
    axis[1].plot(N_array, nu_list, label='Рассчитанное значение')
    axis[1].set_title('График зависимости\nоценок коэффициента вариации от объема выборки')
    axis[1].set_xlabel('Размер выборки')
    axis[1].set_ylabel('Оценка коэффициента вариации')
    axis[1].legend()
    plt.show()


def erlang_model(k, lmbda, N=10000):
    '''Функция моделирования Эрланговского распределения'''
    # Текущая оценка интенсивности и коэффициента вариации потока
    lmbda_old, nu_old = 1, 1
    # Списки оценок интенсивности и коэффициента вариации потока для разных N
    lmbda_list, nu_list = [], []
    # Список значений выборки эрланговского распределения
    u_list = []
    # Счетчик цикла
    loop_cnt = 0
    # Список зничений объема выборки при которых производились
    # расчеты интенсивности и коэффициента вариации потока
    N_array = []

    while True:
        # Добавляем N случайных значений эрланговского закона распределения
        u_list.extend(random_erlang_val(k, lmbda, N * 2**loop_cnt))
        # Добавляем значение нового объема выборки
        N_array.append(len(u_list))

        # Оценка Мат. ожидания и Дисперсии
        me = np.mean(u_list)
        sigma = np.std(u_list)
        # Оценка новых значений интенсивности и коэффициента вариации потока
        lmbda_new = 1/me
        nu_new = sigma/me
        # Добавление новых отценок в списки
        lmbda_list.append(lmbda_new)
        nu_list.append(nu_new)

        # Проверка достаточности выборки
        if abs((lmbda_new - lmbda_old) / lmbda_old) <= 0.01 and abs((nu_new - nu_old) / nu_old) <= 0.01:
           break

        # Перенос новых оценок на место старых
        lmbda_old = lmbda_new
        nu_old = nu_new
        loop_cnt += 1 if len(u_list) != N else 0

    # Построение графика по сформированному распределению
    ploting_erlang_distribution(k, lmbda, data=u_list, add_hist=True)

    # Построение графика зависимости оценок интенсивности и коэффициента вариации потока
    ploting_lambda_nu(k, lmbda, lmbda_list, nu_list, N_array)
    return u_list, lmbda_list, nu_list


def main():
    k, lmbda = 1, 13
    ploting_erlang_distribution(k, lmbda, N=10000)
    erlang_model(k, lmbda, N=1)
    return 0

if __name__ == "__main__":
	main()