Matan2

#!/usr/bin/env python
# coding: utf-8

# # Zadan - 2
# Задание №1. В результате опыта была получена выборочная совокупность.
#
# 1.По данной таблице составить интервальный вариационный ряд, разбив всю вариацию на 8-10 интервалов.
# 2.По сгруппированным данным построить:
#   полигон относительных частот;
#   гистограмму относительных частот;
#   график эмпирической функции распределения.
# 3.Найти числовые характеристики выборки: М_О^*,М_е^*,х ̅_в,σ_(в,) s,V^*,a_s^*,ε_k^*.
# 4. Построить:
#   на чертеже гистограммы её теоретический аналог ƒ(x);
#   на чертеже эмпирической функции F*(x) её теоретический аналог F(x).
# 5. По виду гистограммы и эмпирической функции распределения выборки выдвинуть гипотезу о распределении генеральной совокупности.
# 6. Проверить выполнение правил «трёх сигм».
# 7. Применив критерий согласия Пирсона χ2 с заданным уровнем значимости α=0,01, окончательно принять или отвергнуть выдвинутую гипотезу о распределении генеральной совокупности.
# 8. Построить на одном чертеже:
#   полигон относительных частот p_i^* и кривую распределения p_i;
#   гистограмму теоретических вероятностей (относительных частот) p_i и ƒ(x). Сравнить график p_i с графиком идеально нормального распределения, используя значения  a_s^*,ε_k^*.
# 9. Найти доверительные интервалы для генеральной средней и генерального среднего квадратичного отклонения по уровню надёжности γ=0,99.
#

# # 1. step

# In[2]:


def p(a):
    print(a)
# set arr
s1 = '60 71 62 57 81 55 59 47 75 56 61 60 63 65 59 61 65 58 76 49 '
s2 = '65 64 59 76 58 52 70 77 67 50 65 53 56 64 55 77 51 61 73 64 '
s3 = '45 53 45 58 57 60 48 71 33 65 50 80 58 67 71 51 51 49 66 63 '
s4 = '67 60 67 61 58 36 75 47 68 63 77 75 62 75 70 75 66 53 63 60 '
s5 = '68 67 55 75 71 59 77 58 65 57 55 28 74 71 47 73 40 45 37 66'
str1 = s1 + s2 + s3 + s4 + s5
str1
#str1 = str1.replace(" ", "")

arrInt = []
arrStr = str1.split(' ')
for item in arrStr:
    arrInt.append(int(item.replace(" ", "")))
arrInt_ = arrInt
arrSorted = sorted(arrInt)

print(arrInt)
print(arrSorted)
len(arrInt)


# In[3]:


def setStepH(start, end):
    h =  (end - start) / 10
    x0 = start - h/2
    return h,x0
h, x0 = setStepH(arrSorted[0],arrSorted[-1])

p(x0)
p(h)
p(arrSorted[0])


# In[4]:


diapazon = arrSorted[-1] - x0
p(diapazon)


# In[5]:


h = 6 # меняем шаг для корректировки захвата диаппазона
x0 = arrSorted[0] - h
arrDiap = []
def setArr(x0,h,arrDiap):
    d1 = [x0, x0 + h]
    for i in range(0,10):
        arrDiap.append([(x0+h*i), x0+h*(i+1)])
    return arrDiap
arrDiap = setArr(x0,h,arrDiap)
print('Колличество диапазонов пара чисел -', len(arrDiap))
print('Границы диапазонов -', arrDiap)
print('Минимальное число массива = ',arrSorted[0])
print('Максимальное число массива = ',arrSorted[-1])


# In[6]:


# середина интервала
midlInterval = []
for item in arrDiap:
    midlInterval.append((item[0]+item[1])/2)

print('середина интервала -', midlInterval)


# In[7]:


korteg = list(arrSorted)
#korteg < 37

import numpy as np
import pandas as pd
import math


dfList = pd.DataFrame(arrSorted, columns =['nums'])
# convert numpy array to dataframe
dfList
def getD(dfList, start, end): # вырезает диапазон от и до
    filtered = dfList.loc[(dfList['nums'] <= end) & (dfList['nums'] > start)]
    return list(filtered['nums'].values)
numPak = []
for item in arrDiap:
    numPak.append(getD(dfList, item[0], item[1]))
numPak

chastota = []
for item in numPak:
    chastota.append(len(item))
chastota
print('Частота - ', chastota)
p(sum(chastota))


# In[8]:


relateFreq = []
for i in range(0, len(midlInterval)):
    relateFreq.append(chastota[i]/len(arrSorted))
relateFreq
print('Относительная частота - ', relateFreq)


# In[9]:


idl = range(1,len(relateFreq)+1,1)
T1 = zip(idl,arrDiap, midlInterval, chastota, relateFreq, numPak)
dfT1 = pd.DataFrame(T1, columns =['id','интерв', 'серед интервал','част', 'отн част', 'пакет'])
dfT1.set_index(['id'],inplace = True)

dfT1


# # 2 step

# In[10]:


#!pip install matplotlib
#!pip install pyplot
from matplotlib import pyplot as plt
import numpy as np


# In[ ]:





# In[11]:


PperH = []
for item in relateFreq:
    PperH.append( item/h )
print('Ряд p*/h -', PperH)
XhistList = []
d = 0.00001
for item in arrDiap:
    XhistList.append(item[0]+d)
    XhistList.append(item[1]+d)
XhistList_ = sorted(set(XhistList))
print('Ряд X гистограммы', XhistList_)
#dfT1.hist( column='отн част', bins=8, lw=2,  ec="white" )
#plt.show()
n_bins = XhistList_
x = arrSorted
plt.hist(x, n_bins,lw=2,  ec="white", density=True,histtype='bar', facecolor='b',alpha=0.5)
plt.title('Гистограмма',fontweight ="bold", alpha=1)
plt.show()


# In[12]:


#p(len(XhistList_))
#p(len(relateFreq))
#XhistList_[1:]
plt.plot(XhistList_[1:], relateFreq)
plt.show()


# In[13]:


Fx = []
it = 0
for item in relateFreq:
    Fx.append(it)
    it += item
Fx.append(it)
print('Fx - ', Fx)

plt.plot(XhistList_, Fx)


# In[14]:


print('Fx плотность распределения - ',Fx)


# # 3 step

# 3.1 3.2

# In[58]:


Mmax = max(PperH)
ind = PperH.index(M1)
ML = PperH[ind-1]
MP = PperH[ind+1]
koef = (Mmax - ML)/((Mmax - ML) + (Mmax - MP))
Mo = h*koef + XhistList_[ind]
print('Мода =', Mo)


# In[50]:





MedIAN = []
for i in range(0, len(Fx)-1):
    bool_ = Fx[i]<0.5 and Fx[i+1]>0.5
    if (bool_):
        MedIAN.append(i)
        MedIAN.append(i+1)
# еайдем медиану X и Y
Y = 0.5
k = ( XhistList_[MedIAN[0]] - XhistList_[MedIAN[1]] )  /  ( Fx[MedIAN[0]] - Fx[MedIAN[1]] )
dPlotn = Y - Fx[MedIAN[0]]
dXhistList = dPlotn * k
X =  XhistList_[MedIAN[0]] + dXhistList
print('Координаты медианы = ', X, ' и ', Y)
plt.plot(X, Y, color = "black", marker = "o", linestyle= "--", linewidth = 1, markersize = 10)
plt.plot(XhistList_, Fx)
Moe = [X,Y]
plt.show()


# 3.3

# In[17]:


xp = []
for i in range(0, len(midlInterval)):
    t = midlInterval[i]*relateFreq[i]
    xp.append(t)
print('xp - ', xp)
x2p = []
for i in range(0, len(midlInterval)):
    t = midlInterval[i]**2*relateFreq[i]
    x2p.append(t)
print('x2p - ', x2p)
x3p = []
for i in range(0, len(midlInterval)):
    t = midlInterval[i]**3*relateFreq[i]
    x3p.append(t)
print('x3p - ', x3p)
x4p = []
for i in range(0, len(midlInterval)):
    t = midlInterval[i]**4*relateFreq[i]
    x4p.append(t)
print('x4p - ', x4p)

idl = range(1,len(relateFreq)+1,1)
T1 = zip(idl, midlInterval, chastota, relateFreq, xp,x2p,x3p, x4p)
dfT2 = pd.DataFrame(T1, columns =['id','серед интервал x','част m', 'отн част p','xp','x2p','x3p', 'x4p'])
dfT2.set_index(['id'],inplace = True)

dfT2


# In[18]:


dfT2.sum()


# In[19]:


# выборочное среднее
Xv = dfT2.sum().values[3]
p(Xv)


# In[20]:


# выборочную дисперсию по формуле (8)
Dv = dfT2.sum().values[4] - Xv**2
p(Dv)


# In[21]:


# выборочное среднее квадратическое отклонение (9):
sigmaV = Dv**0.5
p(sigmaV)


# In[22]:


#  исправленную выборочную дисперсию (23):
n = dfT2.sum().values[1]
s2 = n/ (n-1) * Dv
s = s2 ** 0.5
p(s)


# 3.4

# In[23]:


#4) Для характеристики меры колеблемости изучаемого признака относительно выборочной средней вычислим коэффициент вариации по формуле (12):
V = sigmaV / Xv * 100
print(V,' %')


# In[24]:


# Найдём начальные (14) и центральные (15) моменты, используя результаты вспомогательной таблицы:
v1 = dfT2.sum().values[3]
v2 = dfT2.sum().values[4]
v3 = dfT2.sum().values[5]
v4 = dfT2.sum().values[6]
p(v1)
p(v2)
p(v3)
p(v4)


# In[25]:


mu3 = 2*v1**3 - 3*v1*v2 + v3
mu4 = v4 - 4*v1*v3 + 6*v1**2 *v2 - 3*v1**4
p(mu3)
p(mu4)


# In[26]:


#Выборочный коэффициент асимметрии вычислим по формуле (16):
As = mu3 / sigmaV**3
p(As)


# In[27]:


# выборочный коэффициент эксцесса по формуле (17):
ek = mu4 / sigmaV**4 - 3
p(ek)


# # 4 step

# In[28]:


from sympy import *
from sympy.plotting import plot
from sympy.solvers.inequalities import solve_univariate_inequality
init_printing(use_unicode=False, wrap_line=False, no_global=True)
#x = Symbol('x')
a = Xv
sig = s
#stepen = - (x-a)**2 / ( 2*sig**2 )
#top = 1
#bottom = (2*math.pi)**0.5 * sig
#fx1 = top / bottom * exp(stepen)
#aa = dfT1.hist( column='отн част', bins=8, lw=2,  ec="white" )
#plot(fx1, (x, 0, 110))
######1
xxx = range(1,120,1)
yyy = []
def GR(x,a,sig,newArr):
    for item in x:
        stepen = - (item-a)**2 / ( 2*sig**2 )
        top = 1
        bottom = (2*math.pi)**0.5 * sig
        inn = top / bottom * exp(stepen)
        newArr.append(inn)
GR(xxx,a,sig,yyy)
plt.plot(xxx, yyy, color = "red", marker = "o", linestyle= "--", linewidth = 1, markersize = 1)
#######2
n_bins = XhistList_
x = arrSorted
plt.hist(x, n_bins,lw=2,  ec="white", density=True,histtype='bar', facecolor='b',alpha=0.5)
plt.title('Гистограмма',fontweight ="bold", alpha=0.5)
#########2
plt.show()


# In[29]:


def GetLaplas(L):
    bool_ = df.loc[ (df['х']<abs(L)) ]
    return list(bool_['Ф(х)'])[-1]*L/abs(L)
df = pd.read_excel('Individualnaya_rabota__2_primer_varianty_prilozhenia.xlsx', sheet_name='Table 21', header=0)
df # загрузка таблицы лапласа Ф
a = Xv
sig = s
XhistList_
Ylaplas = []
for item in XhistList_:
    itt = ((item) - a)/sig
    out = 0.5 + GetLaplas(itt)
    #print(itt, GetLaplas(itt))
    Ylaplas.append(out)
Ylaplas
print('X нашей функции - ', XhistList_)
print('Fx по Лапласу - ', Ylaplas)
print('Сигма =', sig)
print('a =', a)


# In[30]:


df
plt.plot(XhistList_, Fx)
plt.plot(XhistList_, Ylaplas)
plt.show()


# # 5 step

# По виду гистограммы и функции * F x( ) выдвигаем основную (нулевую) гипотезу H0

# «Генеральная совокупность распределена по нормальному закону с параметрами a = 60.64 , Сигма = 11.08053165774014»

# # 6 step

# Проверим выполнение правила «трёх сигм»:

# In[31]:


Diapazon = [ (a-3*sig), (a+3*sig) ]


# In[32]:


print('НАш диапазон - ', arrSorted)
plt.hist(arrSorted, Diapazon)
print('Вывод все значения попали в диапазон')


# # 7 step

# Проверим соответствие выдвинутой гипотезы H0 опытным данным. Для этого необходимо вычислить теоретические вероятности i
# p и выравнивающие частоты m np i i ′ = .

# In[33]:


#Необходимым условием применения критерия Пирсона является наличие в каждом из интервалов не менее 5 наблюдений
idl = range(1,len(relateFreq)+1,1)
T3 = zip(idl,arrDiap, chastota)
dfT3 = pd.DataFrame(T3, columns =['id','интерв','част'])
dfT3.set_index(['id'],inplace = True)
arrDiap_= []; chastota_ =[];

arrDiap_ = arrDiap[3:]
arrDiap_[0] = [22,46]
chastota_ = chastota[3:]
chastota_[0] = 8
p(arrDiap_)
p(chastota_)


# In[34]:


print('До есть меньше 5')
dfT3


# In[35]:


idl = range(1,len(relateFreq)+1,1)
T3_ = zip(idl,arrDiap_, chastota_)
dfT3_ = pd.DataFrame(T3_, columns =['id','интерв','част'])
dfT3_.set_index(['id'],inplace = True)
p('после все больше 5')
p(arrDiap_)
p(chastota_)
dfT3_


# Искомые вероятности вычисляем по формуле

# In[36]:


def GetLaplas(L):
    bool_ = df.loc[ (df['х']<abs(L)) ]
    #print(bool_)
    return list(bool_['Ф(х)'])[-1]*L/abs(L)
df = pd.read_excel('Individualnaya_rabota__2_primer_varianty_prilozhenia.xlsx', sheet_name='Table 21', header=0)
df # загрузка таблицы лапласа Ф
a = Xv
sig = s

XhistListNew = []
d = 0
for item in arrDiap_:
    XhistListNew.append(item[0]+d)
    XhistListNew.append(item[1]+d)
XhistListNew = sorted(set(XhistListNew))

YlaplasNew = []
for i in range(1,len(XhistListNew)):
    itt1 = ((XhistListNew[i]) - a)/sig
    itt2 = ((XhistListNew[i-1]) - a)/sig
    out = GetLaplas(itt1) - GetLaplas(itt2)
    YlaplasNew.append(out*100)
YlaplasNew
print('X нашей функции - ', XhistListNew)
print('Fx по Лапласу - ', YlaplasNew)
print('Сигма =', sig)
print('a =', a)


# Для дальнейших расчётов заполним вспомогательную таблицу:

# In[37]:


last = []
for i in range(0, len(YlaplasNew)):
    dd =  ( chastota_[i] - YlaplasNew[i] )**2 / YlaplasNew[i]
    last.append(dd)
idl = range(1,len(XhistListNew)+1,1)
T4 = zip(idl, arrDiap_, chastota_, YlaplasNew, last)
dfT4 = pd.DataFrame(T4, columns =['id','интерв','част m','теор част p','(m-m)^2 / m'])
dfT4.set_index(['id'],inplace = True)
p('после все больше 5')
p(YlaplasNew)
p(chastota_)
dfT4


# In[38]:


Hnabl = sum(last)
print('Hnabl = ',Hnabl)
alf = 0.05
interval = len(XhistListNew)
k = interval - 3
Hkr = [alf, k]
print('Hkr - ',Hkr)
print('Согласно таблице приложения 3 критических точек распределения Hkr =', 11.1)


# Т. к. χ2 < χ2кр
# , то нет оснований отвергнуть проверяемую нулевую гипо-
# тезу. Т. е. принимаем предположение, что статистические данные рас-
# пределены по нормальному закону с параметрами a = 86,4 и σ = 8,41.

# In[39]:


teorRad = []
for i in range(0, len(YlaplasNew)):
    round(YlaplasNew[i])*[(arrDiap_[i][0]+0.01)]
    teorRad = teorRad + round(YlaplasNew[i]*100)*[(arrDiap_[i][0]+0.01)]
#print('Возможный ряд для новой теоретической частоты',teorRad)
dfT4.sum()


# По таблице приложения 3 критических точек распределения χ2 , по
# заданному уровню значимости α = 0,05 и числу степеней свободы
# k = 8 − 3 = 5 найдём критическое значение χ2кр(α ;k) = χ2кр(0,05;5) = 11,1.
# Т. к. χ2 χ2кр
#
# набл < , то нет оснований отвергнуть проверяемую нулевую гипо-
# тезу. Т. е. принимаем предположение, что статистические данные рас-
# пределены по нормальному закону с параметрами a = 60.64 и σ = 11.0805.
#

# # 8 step

# In[40]:


YlaplasNew = []
for i in range(1,len(XhistListNew)):
    itt1 = ((XhistListNew[i]) - a)/sig
    itt2 = ((XhistListNew[i-1]) - a)/sig
    out = GetLaplas(itt1) - GetLaplas(itt2)
    YlaplasNew.append(out)
YlaplasNew
plt.plot(XhistList_[1:], relateFreq)
plt.plot(XhistListNew[1:], YlaplasNew)
plt.show()


# In[41]:


plt.plot(xxx, yyy, color = "red", marker = "o", linestyle= "--", linewidth = 1, markersize = 1)
#######2
n_bins = XhistListNew
x = teorRad
plt.hist(x, n_bins,lw=2,  ec="white", density=True,histtype='bar', facecolor='b',alpha=0.5)
plt.title('Стало после перерасчета',fontweight ="bold", alpha=0.8)
#########2
print('Новый ряд', XhistListNew)
plt.show()


# In[42]:


plt.plot(xxx, yyy, color = "red", marker = "o", linestyle= "--", linewidth = 1, markersize = 1)
#######2
n_bins = XhistList_
x = arrSorted
plt.hist(x, n_bins,lw=2,  ec="white", density=True,histtype='bar', facecolor='b',alpha=0.5)
plt.title('Было до теор перерасчёта',fontweight ="bold", alpha=0.9)
#########2
plt.show()


# # 9 step

# ДАНО:<br>
#    -- Применив критерий согласия Пирсона χ2 с заданным уровнем значимости
# α = 0,05, окончательно принять или отвергнуть выдвинутую гипотезу о распре-
# делении генеральной совокупности. <br>
#    -- Найти доверительные интервалы для генеральной средней и генерального
# среднего квадратического отклонения по уровню надёжности γ = 0,99 .

# In[43]:


alf = 0.05
gamma = 0.99
print('Hnabl = ',Hnabl)
alf = 0.05
interval = len(XhistListNew)
k = interval - 3
Hkr = [alf, k]
print('Hkr - ',Hkr)
Hkr_ = 11.1
print('Согласно таблице приложения 3 критических точек распределения Hkr =', Hkr_)


# Значение tγ = t(γ ,n) находим в таблице приложения 6 по заданному
# уровню надёжности

# In[44]:


tgamma = [0.99, 100] # находим по таблице
print('tgamma - ',tgamma)
tgamma_ = 2.627
print('Согласно таблице таблице приложения 6 по заданному уровню надёжности tgamma =', tgamma_)

e = tgamma_ / len(arrInt)**0.5 * sig
print('e =', e)
a = Xv
sig = s


# In[45]:


Xv_m_e = Xv - e
Xv_p_e = Xv + e
DoveritInterval1 = [ Xv_m_e, Xv_p_e ]
print('Доверительный интервал для Хг -', DoveritInterval1)


# In[46]:


q = [0.99, 100] # находим по таблице
print('q - ', q)
q_ = 0.198
print('По таблице приложения 7 величина q(γ ,n) q =', q_)


# In[47]:


sig_mq = sig - sig*q_
sig_pq = sig + sig*q_
# Доверительный интервал для генерального среднего квадратического отклонения
DoveritInterval2 = [ sig_mq, sig_pq ]
print('Доверительный интервал q(γ ,n)  - для Сигма г', DoveritInterval2)