Basic chi2 statistical test with python

chi2.py

import random

class StatisticalTest:
    def __init__(self,
        D_M, D_NM, H_M, H_NM
    ):
        
        self.D_M = D_M
        self.D_NM = D_NM
        self.H_M = H_M
        self.H_NM = H_NM
        
        self.basic_stats()
        self.get_expected()
        self.chi2 = self.get_chi2(
            self.D_M, self.D_NM, self.H_M, self.H_NM,
            self.E_D_M, self.E_D_NM, self.E_H_M, self.E_H_NM
        )
        
    def basic_stats(self,):
        self.TOTAL = self.D_M + self.D_NM + self.H_M + self.H_NM
        self.TOTAL_M = self.D_M + self.H_M
        self.TOTAL_NM = self.D_NM + self.H_NM
        self.TOTAL_D = self.D_M + self.D_NM
        self.TOTAL_H = self.H_M + self.H_NM
        self.p_M = (self.D_M + self.H_M)/self.TOTAL
        
        self.have_mutation = [True] * self.TOTAL_M + [False] * self.TOTAL_NM
        self.have_disease = [True] * self.TOTAL_D + [False] * self.TOTAL_H
        
    def get_expected(self,):
        self.E_D_M = self.TOTAL_D * self.p_M
        self.E_H_M = self.TOTAL_H * self.p_M
        self.E_D_NM = self.TOTAL_D * (1-self.p_M)
        self.E_H_NM = self.TOTAL_H * (1-self.p_M)


    @staticmethod
    def get_chi2(
        O_1, O_2, O_3, O_4,
        E_1, E_2, E_3, E_4,
    ):
        
        diff = (
            (O_1 - E_1)**2 / E_1 +
            (O_2 - E_2)**2 / E_2 +
            (O_3 - E_3)**2 / E_3 +
            (O_4 - E_4)**2 / E_4
        )
        
        return diff
    
    def random_trial(self,):
    
        random.shuffle(self.have_disease)
        random.shuffle(self.have_mutation)

        R_D_M   = sum(k and l for k,l in zip(self.have_disease, self.have_mutation))
        R_D_NM  = sum(k and not l for k,l in zip(self.have_disease, self.have_mutation))
        R_H_M   = sum(not k and l for k,l in zip(self.have_disease, self.have_mutation))
        R_H_NM = sum(not k and not l for k,l in zip(self.have_disease, self.have_mutation))

        chi2_random = self.get_chi2(
            R_D_M, R_D_NM, R_H_M, R_H_NM,
            self.E_D_M, self.E_D_NM, self.E_H_M, self.E_H_NM
        )
        
        return chi2_random >= self.chi2

    def random_trials(self, N):
        
        s = sum(self.random_trial() for _ in range(N))
        
        return s/N

'''        
 D_M = Have Disease with Mutation
 D_NM = Have Disease without Mutation
 H_M = Healthy with Mutation
 H_NM = Healthy without Mutation
 
 Αν τρέξουμε το παρακάτω θα δούμε ότι αν μία μετάλλαξη δεν συσχετίζεται 
 με μία ασθένεια τότε η πιθανότητα να δούμε τη κατανομή:
 D_M=40, D_NM=100, H_M=20, H_NM=100 είναι ~0.028
 
 Αφού το 0.028 < 0.05 λέμε ότι αυτή η κατανομή παραείναι σπάνια 
 για να μην "παίζει κάτι" (δηλαδή μη τυχαίο)
 άρα καταρρίπτουμε την αρχική μας υπόθεση ότι η μετάλλαξη δεν συσχετίζεται
 με την ασθένεια. 
 
 Πειραματιστείτε αλλάζοντας τις τιμές των D_M, D_NM, H_M, H_NM και δείτε
 πως αλλάζει το αποτέλεσμα.
 
 Μπορείτε να σκεφτείτε ότι:
 D_M = Φορές που δεν παίρνω ομπρέλα και βρέχει
 D_NM = Φορές που δεν παίρνω ομπρέλα κα δεν βρέχει
 H_M = Φορές που παίρνω ομπρέλα και βρέχει 
 H_NM = φορές που παίρνω ομπρέλα και δεν βρέχει 
 
 αν το αποτέλεσμα που βγάζει αυτή η ανάλυση < 0.05 τότε ο Θεός με μισεί!
 
 
'''

s = StatisticalTest(
    D_M=40, D_NM=100, H_M=20, H_NM=100
    
)
s.random_trials(10_000)

Το σημαντικό είναι να καταλάβουμε πως με πολυ απλά μαθηματικά (λυκείου) και απλό προγραμματισμό μπορούμε να βγάλουμε κάποιο ασφαλές συμπέρασμα από τα δεδομένα μας.