import random
import math
import numpy as np
import pandas as pd

alpha = 0.05
from google.colab import drive
drive.mount('/content/drive')
file_path = '/content/drive/MyDrive/table.txt'

def PoissonCalculate(lambdaP, x):
        return pow(lambdaP, x) * math.exp(-lambdaP) / math.factorial(x)
def PoissonStandartGenerator(lambdaP):
    r = random.random()
    x = 0
    p = PoissonCalculate(lambdaP, x)
    while(r > p):
        x += 1
        p += PoissonCalculate(lambdaP, x)
    return x

def PoissonModifiedGenerator(lambdaP):
    r = random.random()
    x = 0
    p = math.exp(-lambdaP)
    sumP = p
    while(r > sumP):
        x += 1
        p = p * lambdaP / x
        sumP += p
    return x

xiTable = pd.read_csv(file_path, sep="\t", header=None)
print (xiTable)
xiTable = np.asarray(xiTable, dtype=float)

"ввод:      0         1         2         3         4         5         6         7
0    0   0.95000   0.75000   0.25000   0.10000   0.05000   0.01000   0.00500
1    1   0.00393   0.10153   1.32330   2.70554   3.84146   6.63490   7.87944
2    2   0.10259   0.57536   2.77259   4.60517   5.99146   9.21034  10.59663
3    3   0.35185   1.21253   4.10834   6.25139   7.81473  11.34487  12.83816
4    4   0.71072   1.92256   5.38527   7.77944   9.48773  13.27670  14.86026
5    5   1.14548   2.67460   6.62568   9.23636  11.07050  15.08627  16.74960
6    6   1.63538   3.45460   7.84080  10.64464  12.59159  16.81189  18.54758
7    7   2.16735   4.25485   9.03715  12.01704  14.06714  18.47531  20.27774
8    8   2.73264   5.07064  10.21885  13.36157  15.50731  20.09024  21.95495
9    9   3.32511   5.89883  11.38875  14.68366  16.91898  21.66599  23.58935
10  10   3.94030   6.73720  12.54886  15.98718  18.30704  23.20925  25.18818
11  15   7.26094  11.03654  18.24509  22.30713  24.99579  30.57791  32.80132
12  20  10.85081  15.45177  23.82769  28.41198  31.41043  37.56623  39.99685
13  25  14.61141  19.93934  29.33885  34.38159  37.65248  44.31410  46.92789
14  30  18.49266  24.47761  34.79974  40.25602  43.77297  46.97924  50.89218"

def findTableXi(k, alpha):
    column = 0
    for i in range (1, xiTable.shape[1]):
        if xiTable[0][i] == alpha:
            column = i
    if column == 0:
        return 'Нет нужного значения альфа в таблице'
    i = 1
    while i < len(xiTable) - 1 and xiTable[i][0] < k:
        i += 1
    if k == xiTable[i][0]:
        xi = xiTable[i][column]
    elif k > xiTable[len(xiTable) - 1][0]:
        xi = xiTable[len(xiTable) - 1][column] + (xiTable[len(xiTable) - 1][column] - xiTable[len(xiTable) - 2][column]) / (xiTable[len(xiTable) - 1][0] - xiTable[len(xiTable) - 2][0]) * (k - xiTable[len(xiTable) - 1][0])
    else:
        xi = xiTable[i - 1][column] + (xiTable[i][column] - xiTable[i - 1][column]) / (xiTable[i][0] - xiTable[i - 1][0]) * (k - xiTable[i - 1][0])
    return xi

lambdaP = 0
while lambdaP <= 0:
    lambdaP = float(input('Введите лямбда для распределения Пуассона: '))
    if lambdaP <= 0:
        print('Лямбда не может быть не положительным!')
N = int(input('Введите объём выборки: '))

"ввод: Введите лямбда для распределения Пуассона: 3.75
Введите объём выборки: 1000"

R = np.zeros(N)
for i in range(0, N):
    R[i] = PoissonStandartGenerator(lambdaP)
unique, counts = np.unique(R, return_counts=True)
distribution = np.vstack((unique, counts))
distribution = np.asarray(distribution, dtype=int)
distribution

"ввод: array([[  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11],
       [ 25,  77, 154, 204, 212, 133,  99,  50,  30,  12,   3,   1]])"

R = np.zeros(N)
for i in range(0, N):
    R[i] = PoissonModifiedGenerator(lambdaP)
unique, counts = np.unique(R, return_counts=True)
distribution = np.vstack((unique, counts))
distribution = np.asarray(distribution, dtype=int)
distribution

"ввод: array([[ 60,  66,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,  78,
         79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,  91,
         92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103, 104,
        105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
        118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130,
        131, 132, 133, 134, 137, 138],
       [  1,   1,   3,   3,   4,   2,   3,   9,  14,  20,  19,  32,  30,
         37,  46,  57,  84,  95, 101, 120, 146, 154, 216, 237, 232, 285,
        286, 334, 360, 352, 359, 359, 399, 369, 412, 398, 401, 365, 381,
        355, 332, 297, 271, 265, 250, 216, 201, 175, 145, 139, 102,  83,
         85,  69,  62,  41,  33,  29,  34,  25,  21,   5,  12,   8,   6,
          3,   2,   3,   3,   1,   1]])"

minNumber = 5
expectedDistribution = np.zeros((2, len(unique)))
validArray = True
for i in range (0, len(unique)):
    expectedDistribution[0][i] = distribution[1][i]
    expectedDistribution[1][i] = N * PoissonCalculate(lambdaP, distribution[0][i])
    if expectedDistribution[1][i] < minNumber:
        validArray = False
while not validArray:
    validArray = True
    for i in range (0, expectedDistribution.shape[1] - 1):
        if expectedDistribution[1][i] < minNumber:
            expectedDistribution[0][i + 1] += expectedDistribution[0][i]
            expectedDistribution[1][i + 1] += expectedDistribution[1][i]
            expectedDistribution[0][i] = 0
            expectedDistribution[1][i] = 0
    if expectedDistribution[1][expectedDistribution.shape[1] - 1] < minNumber:
        expectedDistribution[0][expectedDistribution.shape[1] - 2] += expectedDistribution[0][expectedDistribution.shape[1] - 1]
        expectedDistribution[1][expectedDistribution.shape[1] - 2] += expectedDistribution[1][expectedDistribution.shape[1] - 1]
        expectedDistribution[0][expectedDistribution.shape[1] - 1] = 0
        expectedDistribution[1][expectedDistribution.shape[1] - 1] = 0
    expectedDistribution = expectedDistribution[expectedDistribution != 0]
    expectedDistribution = expectedDistribution.reshape(2, int(len(expectedDistribution) / 2))
    for i in range (0, expectedDistribution.shape[1]):
        if expectedDistribution[1][i] < minNumber:
            validArray = False
print(expectedDistribution)


"ввод: [[ 25.          77.         154.         204.         212.
  133.          99.          50.          30.          16.        ]
 [ 23.51774586  88.19154696 165.35915055 206.69893819 193.78025455
  145.33519091  90.83449432  48.66133624  22.81000136  14.28325133]]"

MT = lambdaP
DT = lambdaP
ME = 0
DE = 0
for i in range (0, len(unique)):
    ME += distribution[0][i] * distribution[1][i]
ME /= N
for i in range (0, N):
    DE += (R[i] - ME) * (R[i] - ME)
DE /= (N - 1)
print('Теоретическое значение мат. ожидания: ' + str(MT))
print('Расчётное значение мат. ожидания: ' + str(ME))
print('Теоретическое значение дисперсии: ' + str(DT))
print('Расчётное значение мат. дисперсии: ' + str(DE))

"ввод: Теоретическое значение мат. ожидания: 3.75
Расчётное значение мат. ожидания: 3.843
Теоретическое значение дисперсии: 3.75
Расчётное значение мат. дисперсии: 3.8081591591591484"

intervalsNumber = len(unique)
k = intervalsNumber - 1
XSquaredE = 0
for i in range (0, intervalsNumber):
    XSquaredE += (distribution[1][i] - N * PoissonCalculate(lambdaP, distribution[0][i])) * (distribution[1][i] - N * PoissonCalculate(lambdaP, distribution[0][i])) / (N * PoissonCalculate(lambdaP, distribution[0][i]))
XSquaredT = findTableXi(k, alpha)
print('Теоретическое значение хи-квадрат: ' + str(XSquaredT))
print('Расчётное значение хи-квадрат: ' + str(XSquaredE))

"ввод: Теоретическое значение хи-квадрат: 19.64479
Расчётное значение хи-квадрат: 8.90916222592578"

intervalsNumber = expectedDistribution.shape[1]
k = intervalsNumber - 1
XSquaredE = 0
for i in range (0, intervalsNumber):
    XSquaredE += (expectedDistribution[0][i] - expectedDistribution[1][i]) * (expectedDistribution[0][i] - expectedDistribution[1][i]) / (expectedDistribution[1][i])
XSquaredT = findTableXi(k, alpha)
print('Теоретическое значение хи-квадрат: ' + str(XSquaredT))
print('Расчётное значение хи-квадрат: ' + str(XSquaredE))

"ввод: Теоретическое значение хи-квадрат: 16.91898
Расчётное значение хи-квадрат: 8.332764871150305"