from scipy.stats import poisson
from math import sqrt
from scipy.stats import norm
from scipy.stats import chi2
import numpy as np
from numpy import var
from numpy import mean
from numpy import abs
import random
import time
from random import randint
alpha = 0.05 #hladina testu
lambd = 1  #parametr lambda
poz=0
count = 50 #velikost vzorku nahodneho vyberu
numb =10000  #velikost poctu behu M, MM, W2, CM, KS a FD
cnum = 1000000  #velikost poctu behu testu Qc, Tns a Tan
ano = 0
anocel = 0
pocetb = 99 #pocet testovacich vzorku
n = count
nn = count
cas = 0
def genpomix(a,b,p): #smes Poissonovych distribuci
    aa=[]
    c=[]
    for i in range(count):
        q = random.uniform(0, 1)
        if i==0:
            if q>=p:
                c=poisson.rvs(a,size=1)

                d=c[0]
                d = float(d)
                aa= [d]
            if q<p:
                c=poisson.rvs(b,size=1)

                d = c[0]
                d = float(d)
                aa = [d]
        if i>0:
            if q>=p:
                c = poisson.rvs(a, size=1)
                d = c[0]
                d = float(d)
                aa.append(d)
            if q<p:
                c = poisson.rvs(b, size=1)
                d = c[0]
                d=float(d)
                aa.append(d)
    ma = np.asarray(aa)
    return ma
def genpozero(a,p): #nulou modifikovany Poisson
    aa=[]
    c=[]
    for i in range(count):
        q = random.uniform(0, 1)
        if i==0:
            if q>=p:
                c=poisson.rvs(a,size=1)

                d=c[0]
                d = float(d)
                aa= [d]
            if q<p:
                aa = [0]
        if i>0:
            if q>=p:
                c = poisson.rvs(a, size=1)
                d = c[0]
                d = float(d)
                aa.append(d)
            if q<p:
                aa.append(float(0))
    ma = np.asarray(aa)
    return ma

def genpo(count, lambd):  # generator o velikosti count
    # poisson
    a = poisson.rvs(lambd, size=count)
    #binomicke rozdeleni
    #a = binom.rvs(q, p, size=count)
    #smes Poissonu
    #a = genpomix(q, t, p)
    #nulou modifikovany Poisson
    #a = genpozero(q,p)
    # pro neg. binom.
    #a = nbinom.rvs(q, p, size=count)
    return a


def genpo2(count, c):  # generator o velikosti count, pro bootstrapy
    a = poisson.rvs(c, size=count)
    return a

def test2():  # FD
    anocel = 0
    for i in range(numb):
        ano = 0
        r = genpo(count, lambd)
        prum = r.mean()
        # pro FD statistiku
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
                fd1 = (1 / sqrt(2 * n) * (((n - 1)) * var(r, ddof=1) / mean(r) - n)) * (
                    1 / sqrt(2 * n) * (((n - 1)) * var(r, ddof=1) / mean(r) - n))
        if prum>0:
            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                if prum==0:
                    fd2=0
                else:
                    r2 = genpo2(count, prum)
                    # pro FD statistiku
                    if r2.mean()==0:
                        fd2=0
                    if r2.mean() > 0:
                        fd2 = (1 / sqrt(2 * n) * (((n - 1)) * var(r2, ddof=1) / mean(r2) - n)) * (
                                1 / sqrt(2 * n) * (((n - 1)) * var(r2, ddof=1) / mean(r2) - n))
                if fd1 >= fd2:
                    ano = ano + 1
        if ano > pocetb * (1 - alpha):
            anocel = anocel + 1
    # print(ano)
    print('test zalozeny na FD statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')

def test4():  # w2
    anocel = 0
    for i in range(numb):
        ano = 0
        r = genpo(count, lambd)
        prum = r.mean()
        cm = 0
        d = []
        poz2 = int(prum //2*2 + 5) * 2

        d = []
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            for j in range(poz2):
                d.append(0)
                f = 0
                for k in range(count):
                    if r[k] <= j:  # chybi rovnost
                        f = f + 1
                d[j] = f / count
            for j in range(0, poz2):
                cm = cm + (d[j] - poisson.cdf(j + poz, prum)) * (d[j] - poisson.cdf(j + poz, prum)) * poisson.pmf(j + poz,
                                                                                                                  prum)
            cm1 = count * cm  # nyni cm1 obsahuje hodnotu test. statistiky

            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                r2 = genpo2(count, prum)
                prum2 = r2.mean()
                # pro cm
                bcm = 0
                d2 = []
                for j in range(poz2):
                    d2.append(0)
                    f = 0
                    for k in range(count):
                        if r2[k] <= j:  # chybi rovnost
                            f = f + 1
                    d2[j] = f / count
                j = 0
                for j in range(0, poz2):
                    bcm = bcm + (d2[j] - poisson.cdf(j + poz, prum2)) * (d2[j] - poisson.cdf(j + poz, prum2)) * poisson.pmf(
                        j + poz, prum2)
                cm2 = count * bcm  # nyni cm2 obsahuje hodnotu boostrap. test. statistiky
                if cm1 >= cm2:
                    ano = ano + 1
        if ano > pocetb * (1 - alpha):
            anocel = anocel + 1
    print('test zalozeny na w2 statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')


def test3():  # CM
    anocel = 0
    for i in range(numb):
        ano = 0
        r = genpo(count, lambd)
        prum = r.mean()
        cm = 0
        d = []
        d2 = []
        poz2 = int(prum //2*2 + 5) * 2
        if prum == 0:
            g = randint(1, 20)
            if g == 20:
                anocel = anocel + 1
        else:
            d = []
            for j in range(poz2):
                d.append(0)
                f = 0
                for k in range(count):
                    if r[k] <= j:  # chybi rovnost
                        f = f + 1
                d[j] = f / count
            cm = (d[0] - poisson.cdf(poz, prum)) * (d[0] - poisson.cdf(poz, prum)) * (d[0])
            for j in range(1, poz2):
                cm = cm + ((d[j] - poisson.cdf(poz + j, prum)) * (d[j] - poisson.cdf(poz + j, prum)) * (d[j] - d[j - 1]))
            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                r2 = genpo2(count, prum)
                prum2 = r2.mean()
                d2 = []
                for j in range(poz2):
                    d2.append(0)
                    f = 0
                    for k in range(count):
                        if r2[k] <= j:  # chybi rovnost
                            f = f + 1
                    d2[j] = f / count
                bcm = (d2[0] - poisson.cdf(poz + 0, prum2)) * (d2[0] - poisson.cdf(poz + 0, prum2)) * (
                    d2[0])
                for j in range(1, poz2):
                    bcm = bcm + ((d2[j] - poisson.cdf(poz + j, prum2)) * (d2[j] - poisson.cdf(poz + j, prum2)) * (
                            d2[j] - d2[j - 1]))
                if cm>bcm:
                    ano = ano + 1
        if ano >= pocetb * (1 - alpha):
            anocel = anocel + 1
    print('test zalozeny na CM statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')


def kst(): #KS statistika
    anocel = 0
    for i in range(numb):
        r = genpo(count, lambd)
        prum = r.mean()
        KS1 = 0
        poz2 = int(prum //2*2 + 5) * 2
        F=[]
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            for j in range(poz2):
                F.append(0)
                f = 0
                for k in range(count):
                    if r[k] <= j:  # chybi rovnost
                        f = f + 1
                F[j] = f / count

            for j in range(poz2):
                ks1 = (abs(F[j + poz] - poisson.cdf(j + poz, prum)))
                # print(poisson.cdf(j, prum))
                if ks1 > KS1:
                    KS1 = ks1
            ano = 0
            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                r2 = genpo2(count, prum)
                prum2 = r2.mean()
                KS2 = 0
                F = []
                for j in range(poz2):
                    F.append(0)
                    f = 0
                    for k in range(count):
                        if r2[k] <= j:  # chybi rovnost
                            f = f + 1
                    F[j] = f / count
                for j in range(poz2):
                    ks2 = (abs(F[j + poz] - poisson.cdf(j + poz, prum2)))
                    if ks2 > KS2:
                        KS2 = ks2

                if KS1 >= KS2:
                    ano = ano + 1
            if ano >= pocetb * (1 - alpha):
                anocel = anocel + 1
    print('test zalozeny na KS statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')

def tab(): #pro hodnoty v tabulkach o odhadech m dist. funkce a rozdilu edf
    u=[0,0,0,0,0,0,0,0,0,0]
    u2 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    q = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    for l in range(numb):
        r = genpo(count, lambd)
        prum = r.mean()
        mk = [prum]
        for j in range(poz + 1, poz2 + 1):
            mkc = 0
            for a in range(count):
                mkc = mkc + abs(j - r[a])
            mk.append(mkc / count)
        f0 = (mk[1] + 1 - prum) / 2
        f = [f0]
        F = [f0]
        for i in range(1, poz2):
            f.append((mk[i + 1] - (i + 1 - prum) * (2 * F[i - 1] - 1)) / (2 * (i + 1)))
            F.append(F[i - 1] + f[i])
        for z in range(10):
            u2[z] = u2[z] + F[z]
        for z in range(10):
            x=0
            for w in range(count):
                if r[w] <= z:
                    x = x + 1
            u[z] = u[z] + x
    for z in range(10):
        u2[z]=u2[z]/100

        u[z]=u[z]/5000
        q[z]=u[z]-u[z-1]
    q[0]=u[0]
    #print(u2)
    #print(u)
    #print(q)
    for z in range(10):
        #print(poisson.cdf(z, lambd))
        print(poisson.pmf(z, lambd))



def mtest(): #Mtest
    anocel = 0
    for l in range(numb):
        r = genpo(count, lambd)
        m = 0
        poz = 0

        ano = 0
        prum = r.mean()
        mk = [prum]
        poz2 = int(prum //2*2 + 5) * 2
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            for j in range(poz + 1, poz2 + 1):

                mkc = 0
                for a in range(count):
                    mkc = mkc + abs(j - r[a])
                mk.append(mkc / count)
            # print(mk)

            f0 = (mk[1] + 1 - prum) / 2
            f = [f0]
            F = [f0]
            for i in range(1, poz2):  # od 1
                f.append((mk[i + 1] - (i + 1 - prum) * (2 * F[i - 1] - 1)) / (2 * (i + 1)))
                F.append(F[i - 1] + f[i])

            for j in range(0, poz2):
                F1 = F[j]
                m = m + (F1 - poisson.cdf(poz + j, prum)) * (F1 - poisson.cdf(poz + j, prum)) * (poisson.pmf(j + poz, prum))
            Mn = poz2 * m  # nyni Mn obsahuje hodnotu test. statistiky
            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                r2 = genpo2(count, prum)
                prum2 = r2.mean()
                # pro mtest bootstratz
                mk2 = [prum2]
                for j in range(poz + 1, poz2 + 1):
                    mkc2 = 0
                    for a in range(count):
                        mkc2 = mkc2 + abs(j - r2[a])
                    mk2.append(mkc2 / (count))
                f1 = (mk2[1] + 1 - prum2) / 2
                f2 = [f1]
                F2 = [f1]
                # print(f)
                for i in range(1, poz2):  # od 1
                    f2.append((mk2[i + 1] - (i + 1 - prum2) * (2 * F2[i - 1] - 1)) / (2 * (i + 1)))
                    F2.append(F2[i - 1] + f2[i])
                m1 = 0
                for j in range(poz2):
                    # F2 = Fxk(j, prum2,mk)
                    F22 = F2[j]
                    m1 = m1 + (F22 - poisson.cdf(poz + j, prum2)) * (F22 - poisson.cdf(poz + j, prum2)) * (
                        poisson.pmf(j + poz, prum2))
                Mn2 = poz2 * m1  # nyni cm2 obsahuje hodnotu boostrap. test. statistiky

                if Mn >= Mn2:
                    ano = ano + 1
        if ano > pocetb * (1 - alpha):
            anocel = anocel + 1
    print('test zalozeny na M statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')

def mmtest(): #aka MM test
    anocel = 0
    for l in range(numb):
        r = genpo(count, lambd)
        m = 0
        poz = 0

        ano = 0
        prum = r.mean()
        poz2 = int(prum //2*2 + 5) * 2
        mk = [prum]
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            for j in range(poz + 1, poz2 + 1):

                mkc = 0
                for a in range(count):
                    mkc = mkc + abs(j - r[a])
                mk.append(mkc / count)
            # print(mk)

            f0 = (mk[1] + 1 - prum) / 2
            f = [f0]
            F = [f0]
            for i in range(1, poz2):  # od 1
                #print(F[i - 1], f[i - 1])
                f.append((mk[i + 1] - (i + 1 - prum) * (2 * F[i - 1] - 1)) / (2 * (i + 1)))
                F.append(F[i - 1] + f[i])

            m=(F[0] - poisson.cdf(poz, prum)) * (F[0] - poisson.cdf(poz, prum)) * (F[0])
            for j in range(1, poz2):
                F1 = F[j]
                m = m + (F1 - poisson.cdf(poz + j, prum)) * (F1 - poisson.cdf(poz + j, prum)) * (F[j]-F[j-1])
            Mn = poz2 * m  # nyni Mn obsahuje hodnotu test. statistiky
            for c in range(pocetb):  # 249 bootstrapu Mean dist. t. of Poisson dist. p 255
                r2 = genpo2(count, prum)
                prum2 = r2.mean()
                # pro mtest bootstratz
                mk2 = [prum2]
                for j in range(poz + 1, poz2 + 1):
                    mkc2 = 0
                    for a in range(count):
                        mkc2 = mkc2 + abs(j - r2[a])
                    mk2.append(mkc2 / (count))
                f1 = (mk2[1] + 1 - prum2) / 2
                f2 = [f1]
                F2 = [f1]
                # print(f)
                for i in range(1, poz2):  # od 1
                    f2.append((mk2[i + 1] - (i + 1 - prum2) * (2 * F2[i - 1] - 1)) / (2 * (i + 1)))
                    F2.append(F2[i - 1] + f2[i])

                m1 = (F2[0] - poisson.cdf(poz, prum2)) * (F2[0] - poisson.cdf(poz, prum2)) * (F2[0])
                for j in range(1,poz2):
                    # F2 = Fxk(j, prum2,mk)
                    F22 = F2[j]
                    m1 = m1 + (F22 - poisson.cdf(poz + j, prum2)) * (F22 - poisson.cdf(poz + j, prum2)) * (
                        (F2[j]-F2[j-1]))
                Mn2 = poz2 * m1  # nyni cm2 obsahuje hodnotu boostrap. test. statistiky

                if Mn >= Mn2:
                    ano = ano + 1
            #    print(Mn,Mn2)
        # print(l,ano)
        if ano > pocetb * (1 - alpha):
            anocel = anocel + 1
    print('test zalozeny na MM statistice zamítl nulovou hypotézu v:', anocel, 'z:', numb, 'tj. pro',
          anocel / numb * 100, '%')

def Tantest(): #Zhao Brown, pojmenovano po Anscombem
    anocel = 0
    for i in range(cnum):
        r = genpo(count, lambd)
        y = r *1.0
        prum=r.mean()
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            for j in range(count):
                y[j] = sqrt((r[j] + 3 / 8))
            yprum = mean(y)
            tnew = 0
            for j in range(count):
                tnew = tnew + (y[j] - yprum) * (y[j] - yprum)
            tnew = 4 * tnew
            if tnew <= chi2.ppf(alpha / 2, count - 1) or tnew >= chi2.ppf(1 - alpha/2, count - 1):
                anocel = anocel + 1
    print('testem Tnew neproslo', anocel / cnum * 100, '%')




def Tnstest(): #Neymann Scott
    anocel = 0
    for i in range(cnum):
        r = genpo(count, lambd) * 1.0
        prum = mean(r)
        s = var(r, ddof=1)
        if prum==0:
            g=randint(1,20)
            if g==20:
                anocel=anocel+1
        else:
            if prum>0:
                tns = sqrt((count - 1) / 2) * (s / prum - 1)
            # print(tcc,tc1)

            if tns<norm.ppf(alpha/2) or tns > norm.ppf(1 - alpha/2):
                anocel = anocel + 1
    print('testem Tns neproslo', anocel / cnum * 100, '%')


def Qctest():  # 2. test na str. 277
    anocel = 0
    for i in range(cnum):
        r = genpo(count, lambd) * 1.0
        ex = r.mean()  # prumer
        s = var(r, ddof=1)
        if ex == 0:
            g = randint(1, 20)
            if g == 20:
                anocel = anocel + 1
        else:
            if ex>0:
                q = (count - 1) * s / ex
            if q <= chi2.ppf(alpha / 2, count - 1) or q >= chi2.ppf(1 - alpha/2, count - 1):
                anocel = anocel + 1
    print('testem Test1/Andel neproslo', anocel / cnum * 100, '%')


cas = time.time()

lambd=1
print(lambd, alpha)
#lambd2=lambd #pro smes
#q=lambd    #pro alternativni rozdeleni-musi se jeste upravit funkce genpo co se ma generovat
#t=lambd2
#p #pravdepodobost
#Qctest() #andel
#Tnstest() #NS stat.
#mtest() #Mtest
#mmtest() #MMtest
#test4() # W^2 - cramer von mises
#test3() # CM   - modified cramer von mises#
#kst() # Kolmogorov smirnov
#test2() # Fisher Dispersion

print('konec')
print('doba behu programu:', time.time() - cas)