randvar

import scipy.optimize as opt
from math import erf, pi, sqrt, exp, log
from randtest import sample, expval, disp
import random
from hist import build_hist
import matplotlib.pyplot as plt
import numpy as np

random.seed()


def distrf(x):
    k = x**2/8
    return k * sqrt(2/pi) * exp(-k)


class MaxwellDistribution:
    def __init__(self):
        pass

    def __call__(self, x):
        if x > 0:
            a = 2*sqrt(2)
            # интеграл от плотности
            # при a = 2sqrt(2)
            # (условие нормировки)
            return erf(x/a) - x*exp(-x**2/8)/sqrt(2*pi)
        else:
            return 0.0


def dist_random(dist):
    a = random.random()
    # предполагаем, что вероятность выпадения чисел,
    # больших 20, равна нулю
    # но на самом деле из 100 тыс испытаний
    # самое большое, что я получил --
    # это 11 с копейками
    return opt.bisect(lambda x: dist(x) - a,
                      0.0, 20.0)


if __name__ == "__main__":
    D = MaxwellDistribution()
    N = 10**4

    S = sample(lambda: dist_random(D), N)
    bars = 45
    k = 10 / bars
    # bars = 1 + int(log(N)/log(2))
    G = build_hist(S, bars, 0.0, 10.0)
    M = expval(S)
    V = disp(S, M)

    a = 2 * sqrt(2)

    # Real M,D
    real_M = 4 * sqrt(2/pi)
    real_V = 12 - 32/pi

    print(f' M: {M}, M_отклонение: {abs(real_M-M)/real_M} ')
    print(f' D: {V}, D_отклонение: {abs(real_V-V)/real_V} ')


    plt.style.use('seaborn-dark')
    fig, ax = plt.subplots()

    x = np.fromiter([10*i/bars for i in range(bars)], dtype=float)
    # f = np.vectorize(distrf)
    f = np.vectorize(lambda x: distrf(x)*k)
    f2 = np.vectorize(lambda x: D(x+k)-D(x))
    fs = f(x)
    fr = f2(x)

    x_axis = [f'{X:.1f}' for X in x]
    width = 0.17
    rects1 = ax.bar(x, G, width, label=f'Гистограмма (k={k})')
    plot = ax.plot(x, fs, 'b-', label='kf(x)')
    plot2 = ax.plot(x, fr, 'g-', label='F(x+k)-F(x)')
    ax.set_title('Гистограмма распределения Максвелла')
    ax.set_xticks(x)

    ax.set_ylabel('P')
    ax.set_xticklabels(x_axis)
    plt.xticks(rotation=45)
    ax.legend()
    plt.show()