【C# 单因素方差分析（One Way ANOVA）】

什么是单因素方差分析

单因素方差分析是指对单因素试验结果进行分析，检验因素对试验结果有无显著性影响的方法。
也就是说方差检测是在随机变量Y的不同水平下，检验某个变量X是否有显著性变化。

单因素方差分析是两个样本平均数比较的引伸，它是用来检验多个平均数之间的差异，从而确定因素对试验结果有无显著性影响的一种统计方法。

示例

例如，将抗生素注入人体会产生抗生素与血浆蛋白质结合的现象，以致减少了药效。下表列出了5种常用的抗生素注入到牛的体内时，抗生素与血浆蛋白质结合的百分比。现需要在显著性水平α ＝ 0．05下检验这些百分比的均值有无显著的差异。设各总体服从正态分布，且方差相同。

在这里，试验的指标是抗生素与血浆蛋白质结合的百分比，抗生素为因素，不同的5种抗生素就是这个因素的五个不同的水平。假定除抗生素这一因素外，其余的一切条件都相同。这就是单因素试验。试验的目的是要考察这些抗生素与血浆蛋白质结合的百分比的均值有无显著的差异。即考察抗生素这一因素对这些百分比有无显著影响。这就是一个典型的单因素试验的方差分析问题。

计算步骤

一.计算平均值

二.计算离差平方和

三.计算平均平方

四.方差分析表

F=组间方差/组内方差=MSB/MSE~F(c-1, nT-1)

如果c个总体均值不相等，则组间方差（MSB）会大于组内方差（MSE）。当F值大到某一临界值时，就可以拒绝原假设。临界值的大小由给定的α和自由度决定，一般取α=0.05或者α=0.01。所以，当给定显著性水平为α时，F的拒绝域为F>Fα(c-1,nT-c)。当组间方差和组内方差一样（我们认为这两个组别差别很小），那么F值为1，组间方差远大于组内方差时这个F值也就会比较大。

每一个F值都会对应一个p值，F值越大，p值越小，就越不可能接受原假设，组间差距也就越大。所以，p值越小，拒绝原假设的概率越低，这个特征就越该被保留下来。

C#实现代码

代码实现参考微软官网的讲解，但是在求P值的时候没有提供相应的原函数，在下方列了出来，供参考。
https://docs.microsoft.com/zh-cn/archive/msdn-magazine/2016/october/test-run-anova-with-csharp.

	/*
*  P{ F > F(m, n) } = p
*返回p值
*/
		double sum = 0;//总体数据和
        double SSA=0;//组间离差平方和
        double SSE=0;//组内离差平方和
        double df1 = 0;//组间自由度
        double df2 = 0;//组内自由度
        double f;   //F值
		double p;   //P值
        string F, P;    //用来接收四舍五入后的P、F值
        List<string> Fs = new List<string>();
        List<string> Ps = new List<string>(); 
        int n,m;//二维数组的行列数
        
        private double anova(double[,]arr) {     
            double gmeans;//总均值
            double SST;
            double MSA;
            double MSE;
            n = arr.GetLength(0);   //行数
            m = arr.GetLength(1);   //列数
            double[] total = new double[m];//记录每一组的和
            double[] means = new double[m];//记录每一组的平均值
            df1 = n - 1;//组间自由度 = n-1
            df2 = n * m - n;//n * (m - 1);组内自由度 = n*(m-1)，n为组数，m为每组个数,或者(总数减去组数)
            sum = 0;
            SSE = 0;
            SSA = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    sum = sum + arr[i,j];//求总体和
                }
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    total[i] += arr[i,j];   //求出每组和
                    means[i] = total[i] / m;    //求每组均值
                }
            }

            gmeans = sum / arr.Length;

            for (int i = 0; i < n; i++)
            {
                
                for (int j = 0; j < m; j++)
                {
                    
                    SSE = SSE + Math.Pow((arr[i, j] - means[i]), 2);//组内方差
                }
              
                SSA = SSA + m * Math.Pow((means[i] - gmeans), 2);//组间方差
            }
            SST = SSA + SSE;
            MSA = SSA / df1;//组间均方差（自由度df1 = n-1，n为组数）
            MSE = SSE / df2;//组内均方差（自由度df2 = n*(m-1) 或者 总体数-组数）
            f = MSA / MSE;
            F = f.ToString("f6");   //四舍五入保留6位小数
			p = FDist(f,df1,df2);
            P = p.ToString("f6");   //四舍五入保留6位小数                
            
            return p;
 
        }
        
		/*  下面为用到的一些原函数方法 */
        

        private double FDist(double F, double m, double n)
		{
			double xx;
			double p = 0;

			if (m <= 0 || n <= 0)
			{
				p = -1;
			}
			else if (F > 0)
			{
				xx = F / (F + n / m);
				p = betainc(xx, m / 2, n / 2);
			}
			return (1 - p);

		}

		double betainc(double x, double a, double b)/* 不完全Beta函数 */
		{
			double y, BT, AAA;

			if (x == 0 || x == 1)
				BT = 0;
			else
			{
				AAA = gamma(a + b) - gamma(a) - gamma(b);
				BT = Math.Exp(AAA + a * Math.Log(x) + b * Math.Log(1 - x));
			}
			if (x < (a + 1) / (a + b + 2))
				y = BT * beta_cf(a, b, x) / a;
			else
				y = 1 - BT * beta_cf(b, a, 1 - x) / b;

			return y;
		}

		double beta_cf(double a, double b, double x)
		{
			int count, count_max = 100;
			double eps = 0.0000001;
			double AM = 1;
			double BM = 1;
			double AZ = 1;
			double QAB;
			double QAP;
			double QAM;
			double BZ, EM, TEM, D, AP, BP, AAP, BPP, AOLD;

			QAB = a + b;
			QAP = a + 1;
			QAM = a - 1;
			BZ = 1 - QAB * x / QAP;

			for (count = 1; count <= count_max; count++)
			{
				EM = count;
				TEM = EM + EM;
				D = EM * (b - count) * x / ((QAM + TEM) * (a + TEM));
				AP = AZ + D * AM;
				BP = BZ + D * BM;
				D = -(a + EM) * (QAB + EM) * x / ((a + TEM) * (QAP + TEM));
				AAP = AP + D * AZ;
				BPP = BP + D * BZ;
				AOLD = AZ;
				AM = AP / BPP;
				BM = BP / BPP;
				AZ = AAP / BPP;
				BZ = 1;
				if (Math.Abs(AZ - AOLD) < eps * Math.Abs(AZ)) return (AZ);
			}
			return AZ;
		}

		double gamma(double xx)
		{
			double[] coef_const = new double[7];
			double step = 2.50662827465;
			double HALF = 0.5;
			double ONE = 1;
			double FPF = 5.5;
			double SER, temp, x, y;
			int j;

			coef_const[1] = 76.18009173;
			coef_const[2] = -86.50532033;
			coef_const[3] = 24.01409822;
			coef_const[4] = -1.231739516;
			coef_const[5] = 0.00120858003;
			coef_const[6] = -0.00000536382;

			x = xx - ONE;
			temp = x + FPF;
			temp = (x + HALF) * Math.Log(temp) - temp;
			SER = ONE;
			for (j = 1; j <= 6; j++)
			{
				x = x + ONE;
				SER = SER + coef_const[j] / x;
			}
			y = temp + Math.Log(step * SER);

			return y;
		}