卡方检验值转换为P值

卡方检验作为一种常见的假设检验，在统计学中的地位是显而易见的，如果你还不太清楚可以参看这篇博文:卡方检验用于特征选择，写的非常的浅显易懂，如果你还想再扩展点卡方检验方面的知识，可以参看这篇博文卡方检验基础，写的也很有意思。前辈的功底都很深厚，小弟就就不再阐述卡方检验的原理、意义及如何计算了，理解了其实很简单就那么个公式，再根据实际业务场景关键看你选择哪一个。从chi-squared value 到p-value，相信大多数同学和我一样，查表，因为大学课本上就是这么写的。假如在实际业务场景中，自由度和显著性水准都不确定的情况下，怎么办呢？查表就显得不那么地道了。

这时可能很多同学想到了著名的fisher精确检验，因为这个检验能直接求出的精确的p-value，但是在检验数据样本比较大的情况下，fisher精确检验的计算复杂度会让你显得那么的力不从心，本系列的后面将会讲到fisher精确检验的原理并给出其源码及与chi-squared的效率对比。还是抓紧时间侃侃怎么通过chi-squared计算p-value吧，估计心急的同学就等不及了。ok，咱们攻城师还是用代码说话，先上代码：

public double chisqr2pValue(int dof, double chi_squared) {if (chi_squared < 0 || dof < 1) {return 0.0;}double k = ((double) dof) * 0.5;double v = chi_squared * 0.5;if (dof == 2) {return Math.exp(-1.0 * v);}double incompleteGamma = IncompleteGamma.log_igf(k,v);// 如果过小或者非数值或者无穷if (Math.exp(incompleteGamma) <= 1e-8|| Double.isNaN(Math.exp(incompleteGamma))|| Double.isInfinite(Math.exp(incompleteGamma))) {return 1e-14;}double gamma = Math.log(Gamma.getApproxGamma(k));incompleteGamma -= gamma;if(Math.exp(incompleteGamma) > 1){return 1e-14;}double pValue = 1.0 - Math.exp(incompleteGamma);return (double) pValue;
}

这个chisqr2pValue函数引用到了两个函数，一个为getApproxGamma(k) 一个为log_igf(k,v)，如果你对该函数的实现原理不太清楚，wiki中侃的很清楚：卡方分布，本文也不再阐述卡方分布的密度函数、自由度等等；具体求P-value说白了就是计算卡方分布的分布函数值，如下的公式：

其中，分子为不完全伽马函数，分母为伽马函数；那么上述chisqu2pvalue就是实现了上述公式。Ok,现在的问题转换为怎么求分子与分母。

在两年前我曾为了实现伽马函数的功能代码敲个积分，但是效果不太理想。但，今天发现个利器：斯特灵公式，估计大多数同学跟我一样，第一感觉斯特灵公式不是用来求阶乘的吗？不错，确实是，在实现fisher时我也用到了斯特灵公式；可是哥们确实有点孤陋寡闻，斯特灵还有个求伽马函数的近似公式，如下：

如果你想了解详细的斯特灵公式的推导信息，也是参考wiki:斯特林公式是不是觉得维基太牛逼了？好吧，还是赶紧上代码：

public static double getApproxGamma(double n) {// RECIP_E = (E^-1) = (1.0 / E)double RECIP_E = 0.36787944117144232159552377016147;// TWOPI = 2.0 * PIdouble TWOPI = 6.283185307179586476925286766559;double d = 1.0 / (10.0 * n);d = 1.0 / ((12* n) - d);d = (d + n) *RECIP_E;d = Math.pow(d,n);d *= Math.sqrt(TWOPI/ n);return d;}

好，剩下的就差不完全伽马函数没有解了，怎么求？问大神wiki:不完全伽马函数哎呦，在里面又发现个牛逼的不用计算积分的公式：

而M函数为什么呢？就是传闻中的合连几何函数，如果对这个函数感兴趣的可以继续参考wiki:合连几何函数 ,不过这里我们依旧只使用它的一个公式：

Ok,公式知道了，代码实现就很简单了，只不过在这里为了便于大数的计算，我采用了计算其log值，代码如下：

 public static double log_igf(double s, double z) {if (z < 0.0) {return 0.0;}double sc = (Math.log(z) * s) - z - Math.log(s);double k = KM(s, z);return Math.log(k) + sc;}private static double KM(double s, double z) {double sum = 1.0;double nom = 1.0;double denom = 1.0;double log_nom = Math.log(nom);double log_denom = Math.log(denom);double log_s = Math.log(s);double log_z = Math.log(z);for (int i = 0; i < 1000; ++i) {log_nom += log_z;s++;log_s = Math.log(s);log_denom += log_s;double log_sum = log_nom - log_denom;sum += Math.exp(log_sum);}return sum;}

Ok,简单吧？！以后你再也不用坑爹地查表了。下一篇再介绍下fisher精确检验吧。

卡方检验值转换为P值相关推荐

将RGB值转换为灰度值的简单算法(转)
将RGB值转换为灰度值的简单算法原文:http://blog.163.com/zhaowei0425@126/blog/static/475860302011311103956748/ RGB是如何 ...
matlab rad2deg,在R编程中将弧度值转换为度值–rad2deg（）函数
R语言中的rad2deg()函数用于将指定的弧度值转换为度数值.注意:此函数需要安装"grid"包.语法:rad2deg(x)参数:x:要转换的弧度值示例1:filter_none ...
2021.04.22【RNA-seq流程】丨count值转换为FPKM值优化2.0
优化内容解决每次转换需要设置样本数和基因数目实现基因count值与length精准匹配摘要大概半年前,我写过一篇将HTseq生成的基因COUNT值转换为FPKM值文章,用于对count的入门级 ...
2020.08.14【RNA-Seq流程】丨将HTseq生成的基因COUNT值转换为FPKM值
2021.04.22更新:2.0版本 count转换FPKM介绍公式介绍 R语言脚本公式介绍下方公式是FPKM的综述文章写的,参数并不好理解,但是对数量级表示很明白 https://harold ...
ADC值转换为电压值（机械语言得出电压值）
如何利用单片机的ADC模块(或者独立的ADC芯片)得到接入ADC管脚上的实际电压值? 这个问题,是第一次接触ADC时候,大家都会遇到的问题. 会读到什么值单片机会读到什么值?需要看一个特性,就是几位 ...
C语言之将弧度值转换为角度值
1 #include<stdio.h>2 #include<stdlib.h>3 4 #define pi 3.1415925 6 void conver_radian(flo ...
如何在Microsoft Excel中将文本转换为日期值
Analysis of business data often requires working with date values in Excel to answer questions such ...
不再迷惑，无值和 NULL 值
在关系型数据库的世界中,无值和NULL值的区别是什么?一直被这个问题困扰着,甚至在写TSQL脚本时,心有戚戚焉,害怕因为自己的一知半解,挖了坑,贻害后来人,于是,本着上下求索,不达通幽不罢休的决心(开 ...
【C++grammar】左值、右值和将亡值
目录 C++03的左值和右值 C++11的左值和右值将亡值在C++03中就有相关的概念 C++03的左值和右值通俗的理解: (1) 能放在等号左边的是lvalue (2) 只能放在等号右边的是r ...
医学图像的CT值与像素值总结及转换代码
目录一.CT图像的调窗 1.Window width 2.Window level/center 二.DICOM文件中窗宽窗位对应字段三.CT值与像素值转换(线性映射) 1.itk-snap软件和 ...

卡方检验值转换为P值

卡方检验值转换为P值相关推荐

最新文章

热门文章