import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import datetime
start = datetime.datetime.now()
plt.rcParams['axes.grid'] = True
XX=[None, 11.3, 14.8, 7.6, 10.5, 12.7, 3.9, 11.2, 5.4, 8.5, 5.0, 4.4, 7.3, 2.9, 5.7, 6.2, 7.3, 3.3, 4.2, 5.5, 4.2]
I = [None,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
N=20
    # Задаем G-ю функцию которая включает в себя произведения плотностей вроятности
    # каждого независимого результата Xi в которых МО и СКО изменяются по линейной зависимости
def G(M1, Mn, S1, Sn):
    def loc (M1, Mn, i, n):
        return (M1*(n-i)/(n-1)) + (Mn*(i-1)/(n-1)) 
    def scale (S1, Sn, i, n):
        return (S1*(n-i)/(n-1)) + (Sn*(i-1)/(n-1))
    def fnorm (x):
        return (np.exp((-x**2)/2))/(np.sqrt(2*np.pi))
    def y (M1, Mn, S1, Sn, i, n, x):
        return (x-loc(M1,Mn,i,n))/scale(S1,Sn,i,n)
    def F (M1, Mn, S1, Sn, i, n, x):
        return (fnorm(y(M1, Mn, S1, Sn, i, n, x)))/scale(S1, Sn, i, n)
    # Распаковка значений x и i из глобальных переменных
    x = XX[1:]
    i = I[1:]
    n=N
    # Вычисляем значение функции F для каждого x и i
    values=[F(M1, Mn, S1, Sn, i_val, n, x_val) for i_val, x_val in zip(i, x)]
    
    # Вычисляем произведение всех значений
    result = np.prod(values)
    
    return result
    # находим сомножитель К для получения общей ПВ оценок
options={'epsrel':1e-20, 'epsabs':1e-20}
K, errorK = integrate.nquad(G, ranges=[[0, 30],[0, 10],[0, 10],[0, 10]], opts=[options, options, options, options])
# K = 2.9613332457351404e-18    errorK = 9.999171231431291e-21
    #  формируем ПВ оценок
def pdf(Mn, S1, Sn,       M1):
    return (G(M1, Mn, S1, Sn) / K)
    #  строим график автономной ПВ оценок для параметра М1 (уменьшаем значения ошибок для оперативности)
def pdf_m1 (M1):
    return [(integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(m,), epsabs=0.1, epsrel=0.1)[0]) for m in M1]
x = np.arange(4, 16, 0.2)
plt.plot(x, pdf_m1(x), 'k', lw=3)
plt.ylim(bottom=0)
plt.title('Fm1(m1)')
plt.show()
    # находим несмещенную оценку М1 и её ско sm1 (примерные значения по методу ММП М1 = 9.66)
def F1 (M1):
    return integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(M1,), epsabs=0.01, epsrel=0.01)[0]
M1, _ = integrate.quad(lambda M1: M1 * F1(M1), 4, 16)
print(M1)
Dm1, _ = integrate.quad (lambda x: ((x-M1)**2) * F1(x), 4, 16)
sm1 = np.sqrt(Dm1)
print(sm1)
print("Время вычисления М1 и СКОм1:", datetime.datetime.now()-start)