卡方檢驗（python代碼實現）

python金融風控評分卡模型和數據分析微專業課（博主親自錄制視頻）：http://dwz.date/b9vv

醫葯統計項目QQ：231469242

 

分類變量檢驗方法

 

 

 

 

卡方分布繪圖

 

如果多個符合正態分布的獨立隨機變量z1,z2,z3.....zk，
z1+z2+z3+....z_k呈現卡方分布，自由度k.
有幾個正態分布相加，就有幾個自由度

 

 

# -*- coding: utf-8 -*-
# Toby QQ：231469242
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
import seaborn as sns
import math,pylab,matplotlib,numpy
from matplotlib.font_manager import FontProperties

#設置中文字體
font=FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=15)

n=10

#繪制自由度為n的卡方分布圖,n表示生成卡方數組的個數
def Get_chisquareDatas(n):
    #標准正太分布
    normalDistribution=stats.norm(0,1)
    list_data=[]
    for i in range(n):
        normal_data=normalDistribution.rvs(30)
        chisquare_data=normal_data**2
        list_data.append(chisquare_data)
    return list_data
    
def Plot_chisquare(n):
    list_data=Get_chisquareDatas(n)
    sum_data=sum(list_data)
    plt.hist(sum_data)
    
Plot_chisquare(2)
Plot_chisquare(3)
Plot_chisquare(10)

 

官方繪圖代碼

# -*- coding: utf-8 -*-

from scipy.stats import chi2
import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(1, 1)

df = 20
mean, var, skew, kurt = chi2.stats(df, moments='mvsk')

#繪制函數的起始點和終止點
#pdf為概率密度函數
#百分比函數(PPF) :the inverse of the CDF. PPF  函數和連續分布函數CDF相逆，
#比如輸入哪一個點，可以得到低於等於20的概率？
#ppf(0.01, df)表示輸入哪個點，得到概率低於0.01
initial=chi2.ppf(0.01, df)
end=chi2.ppf(0.99, df)
x = np.linspace(initial,end, 100)

#概率密度函數用於繪圖
ax.plot(x, chi2.pdf(x, df), 'r-', lw=5, alpha=0.6, label='chi2 pdf')
plt.title("df is %d"%df)
plt.show()

 

 

卡方檢驗代碼

可汗學院的問題

# -*- coding: utf-8 -*-
'''
卡方公式(o-e)^2 / e
期望值和收集到數據不能低於5，o(observed)觀察到的數據，e（expected）表示期望的數據
(o-e)平方，最后除以期望的數據e
'''

import numpy as np
from scipy import stats
from scipy.stats import chisquare         
list_observe=[30,14,34,45,57,20]
list_expect=[20,20,30,40,60,30]


std=np.std(data,ddof=1)
print(chisquare(f_obs=list_observe, f_exp=list_expect))
p=chisquare(f_obs=list_observe, f_exp=list_expect)[1]
'''
返回NAN，無窮小
'''

if p>0.05 or p=="nan":
   print"H0 win,there is no difference"
else:
   print"H1 win,there is difference"

 

 

contigency table聯立表

 

 

 

 

 

測試數據

第一行：草本葯1，草本葯2，安慰劑

第二行：生病和非生病

H0：草本葯和疾病無關系

H1：草本葯和疾病有關系

可汗學院計算出來的卡方值2.53；自由度2，顯著性0.1，的關鍵值4.6

卡方值2.53<關鍵值4.6,  H0成立，兩者無關系

 

 

 

python代碼與可汗學院算出結果一致，此版本體現算法推導過程。缺點就是要自己計算出期望值列表

 

# -*- coding: utf-8 -*-
'''
卡方公式(o-e)^2 / e
期望值和收集到數據不能低於5，o(observed)觀察到的數據，e（expected）表示期望的數據
(o-e)平方，最后除以期望的數據e
聯立表contigency table計算
'''

from scipy.stats import chisquare   

list_observe=[34,38,28,50]
list_expect=[29.76,42.2,32.24,45.77]

row=2
colume=2



def Contigency_table(row,colume,list_observe,list_expect):
    degreeFreedom=(row-1)*(colume-1)
    print(chisquare(f_obs=list_observe, f_exp=list_expect,ddof=degreeFreedom))
    p=chisquare(f_obs=list_observe, f_exp=list_expect)[1]


    if p>0.05 or p=="nan":
       print"H0 win,there is no difference"
    else:
       print"H1 win,there is difference"

Contigency_table(row,colume,list_observe,list_expect)

 

此版本不用算出期望值，更加方便，參考的是2*2聯立表，自由度=1，critical value=2.7

# -*- coding: utf-8 -*-

#獨立性檢驗test for independence,也是卡方檢驗chi_square
#前提條件：a,b,c,d 必須大於5

#2.706是判斷標准（90概率），值越大，越有關，值越小，越無關
def value_independence(a,b,c,d):
    if a>=5 and b>=5 and c>=5 and d>=5:
        return ((a+b+c+d)*(a*d-b*c)**2)/float((a+b)*(c+d)*(a+c)*(b+d))

#返回True表示有關
#返回False表示無關
def judge_independence(num_independence):
    if num_independence>2.706:
        print ("there is relationship")
        return True
    else:
        print("there is no relationship")
        return False

a=34
b=38
c=28
d=50
chi_square=value_independence(a,b,c,d)
relation=judge_independence(chi_square)

 

python官網版本，更加方便和科學

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html

import scipy.stats as stats

data = np.array([[43,9],
[44,4]])
V, p, dof, expected = stats.chi2_contingency(data)
print(p)

 

 

 

 

 

python機器學習生物信息學系列課（博主錄制）：http://dwz.date/b9vw