Whole Mat Stat work (Normal Distribution Generator, Confidence Intervals, Student-t, Chi-square Normality Test)

SageMath

import numpy as np
import scipy as sp
import scipy.stats
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import math

@interact
def _(mean1=0, sig1=2, mean2=-1, sig2=2, n=10):
    X1 = np.zeros(n)
    X2 = np.zeros(n)
    for i in range(0, n):
        X1[i] = np.random.normal(mean1, sig1)
        X2[i] = np.random.normal(mean2, sig2)

    #print table([(i, X1[i], X2[i]) for i in [0..n-1]], header_row=['i','X1','X2'], frame=True)            #Whole distribution, both
#==============================================================================================
    #Confidence Interval
    def atv(data):
        buf = 0
        Aver = 0                                              #Average through Volume(n)
        for i in range(0, n):
            buf += data[i]
        return buf/n

    def ci(data, confidence):
        m, se = atv(data), np.std(data)/np.sqrt(n)
        if (confidence == 0.95):
            conf = 1.96
        elif (confidence == 0.90):
            conf = 1.65
        elif (confidence == 0.99):
            conf = 2.57
        h = se * conf                                         #Za/2 * σ/√(n)
        return m, m-h, m+h
    #print 'CI 1: Average = %.3f Minimal = %.3f Maximal = %.3f' % ci(X1, 0.95)
    #print 'CI 2: Average = %.3f Minimal = %.3f Maximal = %.3f' % ci(X2, 0.95)
    #END Confidence Interval
#============================================================================
    #print np.mean(X1)
    n1 = n
    n2 = n
    r = n1 + n2 - 2
    #print 'r =', r
    #3 methods of finding t
    t1 = np.absolute((atv(X1) - atv(X2)))/np.sqrt((n1*(np.var(X1)) + (n2*(np.var(X2))))) * np.sqrt(n1 * n2 / (n1 + n2) * r)
    t2 = sp.stats.ttest_ind(X1, X2)
    t3 = (atv(X1) - atv(X2)) / np.sqrt(sp.stats.sem(X1)**2 + sp.stats.sem(X2)**2)
    #END Student's t-test

    #print 't1 = %.10f' % t1 
    #print 't2 = %.10f p_value = %f' % t2
    #print 't3 = %.10f' % t3

    p_val = sp.stats.ttest_ind(X1, X2)[1]
    #print 'p_val =', p_val

    NoInteresting, Interesting = sp.stats.t.ppf([1-0.025, 1-0.005], r)[0], sp.stats.t.ppf([1-0.025, 1-0.005], r)[1]
    #print NoInteresting, '....', Interesting

    #if (p_val <= NoInteresting):
        #print 'H0 - accepted'
    #elif (NoInteresting < p_val = Interesting):
        #print 'H0 - rejected, H1 - approved'
    #END Student
#===========================================================================
    #Test Normality chi square    
    def chi_square(Obs, Exp):        
        return (np.absolute(Obs - Exp))**2/Exp
    #print chi_square(90, 100) + chi_square(60, 50)
    def chi_sq(o, e):
        return ((o - e)**2)/e

    #Sturgess formula 
    arra = X1                                                                #input what exact array
    Sturgess_delta = max(arra) - min(arra)
    Sturgess_h = Sturgess_delta / (1 + 3.322 * ln(n))
    Sturgess_M = math.floor((Sturgess_delta + Sturgess_h) / Sturgess_h)
    Sturgess_H = Sturgess_delta / (Sturgess_M - 1)
    #END Sturgess formula
    print 'max', max(arra)
    print 'min', min(arra)
    #print 'delta', Sturgess_delta
    #print 'h', Sturgess_h
    print 'M', Sturgess_M
    print 'H', Sturgess_H
    #print np.histogram(X1.astype(int), bins=13)[0]
    #print np.histogram(X1.astype(int), bins=13)[1]

    def intervals(arra):
        global interv
        interv = np.zeros(int(Sturgess_M) + 1)
        interv[0] = min(arra) - (Sturgess_H/2)
        interv[1] = min(arra) + (Sturgess_H/2)
        for i in range(2, int(Sturgess_M) + 1):
            interv[i] = interv[i - 1] + (Sturgess_H)
        return interv
    #print intervals(X1)
    interval = [i for i in intervals(X1)]
    count = np.zeros(Sturgess_M + 1)
    for Ex in X1:                                                       #Zapolnyau viborku v count v opredelennom intervale
        #while Ex > interval in interv:
        for i in range(1, len(intervals(X1))):
            if intervals(X1)[i - 1] < Ex < intervals(X1)[i]:
                #print '\nintervals', intervals(X1)[i - 1]
                #print 'Ex', Ex
                #print 'intervals', intervals(X1)[i]
                count[i - 1] += 1
                break
    #print 'count', count

    #LAPLACCIAN FUNC
    def lapl(x):
        return 2.7182**(-((x**2) / 2))

    def errfunc(t):
        integral = sp.integrate.quad(lapl, 0, t)[0]
        return 1 / np.sqrt(2 * 3.14) * integral
    #END LAPLACCIAN FUNC

    def atd(arr):                                                       #Ispravl vib srednee kvadr otklon
        Sd = np.sqrt(n / (n-1) * np.std(arr))
        return Sd
    def forLapl(arr):
        #for i in intervals(arr):                                      #LOOK HERE!
            #X = (i - atv(arr)) / atd(arr)
        X = [(i - atv(arr)) / atd(arr) for i in intervals(arr)]        #Znacheniya dlya podscheta laplac funcsii
        X[0], X[0-1] = 0, 0
        return X
    #print 'forLalp', forLapl(X1)

    LaplArray = [errfunc(i) for i in forLapl(X1)]
    LaplArray[0], LaplArray[0-1] = -0.5, 0.5
    #print 'LaplArray', LaplArray

    Pi = np.zeros(Sturgess_M)
    for i in range(1, len(LaplArray)):
        Pi[i-1] = LaplArray[i] - LaplArray[i-1]
    #print 'Pi', Pi
    Chastota = [Prob * n for Prob in Pi]
    #print 'Chastota', Chastota
    #print len(interval), len(LaplArray), len(Pi), len(Chastota)
    print table(rows=[([i+1], '%.2f' % interval[i], '%.2f' % interval[i+1], '%.3f' % forLapl(X1)[i], '%.3f' % forLapl(X1)[i+1], '%.3f' % LaplArray[i], '%.3f' % LaplArray[i+1], '%.3f' % Pi[i], '%.3f' %  Chastota[i]) for i in [0..int(Sturgess_M)-1]], header_row=['N', 'Interval[i]', 'Interval[i+1]','forLapl [i]', 'forLapl [i+1]', 'F(forLapl [i])', 'F(forLapl [i+1])', 'Pi', 'Ni* = n*Pi'], frame=True, align='center')
    print 'Summa Ni*', sum(Chastota)

    print table(columns=[('%.3f' % Chastota[i], count[i]) for i in [0..int(Sturgess_M)-1]], header_column=['Ni*', 'Ni'], frame=True, align='center')
    print 'X(experemental)', sum(chi_sq(count[:len(Chastota)], Chastota)), 'Pri %d' % (Sturgess_M-3), 'stepeney svobodi'

    #Xa = atv(X1)
    #Xd = np.sqrt(np.var(X1))
    #Xt = np.zeros(n)
    #for i in range(0, len(Xt)):
    #    Xt[i] = np.random.normal(Xa, Xd)

    #bins = np.arange(0, 10, 1)
    bins = Sturgess_M
    Total = 0

    BuffArr = X1                                         # BUFFER


    #Reshape incomming array to N and then ... NO NEED TO UNCOMMENT
    RS = np.reshape(BuffArr, (n/10,10))
    for SH in range(0, n/10):
        #print RS[SH]
        #print sturgess(X1)[2]
        his = np.histogram(RS[SH].astype(int), bins=bins)[0]
        #print his

        def gaussian(x, a, m, s):
            return a*2.71828^(-(x-m)^2/(2*s^2))
        #print sp.integrate.quad(gaussian, 1.0357, 1.1071, args=(a, m, s))[0]

        #a = 0.266*10
        #m = np.mean(BuffArr)
        #s = np.std(BuffArr)
        #N = 8
        #I = np.zeros(N)

        #for i in range(0, N):
        #    I[i] = sp.integrate.quad(gaussian, i, i+1, args=(a, m, s))[0]

        #Chi = np.zeros(N)
        #for i in range(0, N):
        #    Chi[i] = chi_square(his[i], I[i])
        #SChi = np.sum(Chi)

        Chi_10 = 14.684
        Chi_90 = 4.16
        Chi_75 = 5.89883
        Chi_25 = 11.38875

        #print table(rows=[([i, i+1], '%.2f' % I[i], his[i], '%.3f' % Chi[i]) for i in [0..N-1]], header_row=['Bin','Expected','Obsserved', 'Chi Square test'], frame=True, align='center') 

        #if (SChi = Chi_10):
        #    print '\nH0 - rejected, H1 - approved (Difference permitted)\n'
        #else:
        #    print '\nRetry, not clear\n'

        #elif (Chi_90 < Chi < Chi_10):
        #    print 'Retry, not clear'

        #Total += SChi
    #print Total
#END

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s