统计学：卡方分布临界值的计算

前言

做应用统计分析的课后习题时，有一道题的 $\chi^2$ 分布临界值在表上查不到。老师给出的解决办法是在EXCEL中用=CHISQ.INV.RT( $\alpha$ , n)来计算。但是这样未免有些无聊，为什么不自己写一个程序呢？

准备知识

卡方分布的分布函数

$F(x,k)=\frac{\gamma (\frac{k}{2},\frac{x}{2})}{\Gamma(\frac{k}{2})}$

其中：

$\gamma (s,x)=\int_0^xe^{-t}t^{s-1}dt$

$\Gamma (s)=\int_0^\infty e^{-t}t^{s-1}dt$

临界值

$P(X<x)=F(x,k)$

two tailed z critical value

所谓临界值就是求得一个 $x$ 使得 $P(X<x)$ 大于等于置信水平。
由此我们列出关于 $x$ 的方程，其中 $\alpha$ 为拒绝域面积，1- $\alpha$ 为置信水平：

$F(x,k)=1-\alpha$

scipy模块

SciPy 包含的模块有最优化、线性代数、积分、插值、特殊函数、快速傅里叶变换、信号处理和图像处理、常微分方程求解和其他科学与工程中常用的计算。
本次用到SciPy中的积分模块和求解非线性方程组的模块

含参数定积分：

from scipy import integrate
def f(x, a, b):
    return a * x + b
v, err = integrate.quad(f, 1, 2, args = (-1, 1))
print v

非线性方程求解（最小二乘法）：

1 2	r = fsolve(func,[]) #[]中为初值

python实现

from scipy import integrate
import math
from scipy import optimize
#定义伽马函数
def gamma(s):
    def f(t,s):
        return math.exp(-t)*pow(t,s-1)
    v,err=integrate.quad(f,0,float('inf'),args=(s))
    return v
#定义上不完全伽马函数
def LI_gamma(s,x):
    def f(t,s):
        return math.exp(-t)*pow(t,s-1)
    v,err=integrate.quad(f,0,x,args=(s))
    return v
#定义卡方分布的分布函数
def Fchi(x,k):
    F=LI_gamma(k/2,x/2)/gamma(k/2)
    return F
#方程求解
def cv_chi(alpha,n):
    def func(x):
        return Fchi(x,n)-1+alpha
    r=optimize.fsolve(func,0)
    return r