2

我如何更好地计算定积分？我正在使用一个函数来集成，另一个函数来递归地找到阶乘。

我想改进算法或效率，甚至在这方面的准确性。

    public static double testStatistic(double meanTreatmentSumOfSquares, double meanErrorSumOfSquares)
    {
        return (meanTreatmentSumOfSquares / meanErrorSumOfSquares);
    }

    public static double pValue(double fStatistic, int degreeNum, int degreeDenom)
    {
        double pValue = 0;
        pValue = integrate(0, fStatistic, degreeNum, degreeDenom);

        return pValue;

    }

    public static double integrate(double start, double end, int degreeFreedomT, int degreeFreedomE)
    {
        int iterations = 100000;
        double x, dist, sum = 0, sumT = 0;
        dist = (end - start) / iterations;
        for (int i = 1; i <= iterations; i++)
        {
            x = start + i * dist;
            sumT += integralFunction(x - dist / 2, degreeFreedomT, degreeFreedomE);
            if (i < iterations)
            {
                sum += integralFunction(x, degreeFreedomT, degreeFreedomE);
            }
        }
        sum = (dist / 6) * (integralFunction(start, degreeFreedomT, degreeFreedomE) + integralFunction(end, degreeFreedomT, degreeFreedomE) + 2 * sum + 4 * sumT);
        return sum;
    }

    public static double integralFunction(double x, int degreeFreedomT, int degreeFreedomE)
    {
        double temp=0;
        temp = ((Math.Pow(degreeFreedomE, degreeFreedomE / 2) * Math.Pow(degreeFreedomT, degreeFreedomT / 2)) / (factorial(degreeFreedomE / 2 - 1) * factorial(degreeFreedomT / 2 - 1))) * (factorial(((degreeFreedomT + degreeFreedomE) / 2 - 1)))*((Math.Pow(x, degreeFreedomE / 2 - 1)) / (Math.Pow((degreeFreedomT + degreeFreedomE * x), ((degreeFreedomE + degreeFreedomT) / 2))));
        return temp;
    }

    public static double factorial(double n)
    {
        if (n == 0)
        {
            return 1.0;
        }
        else
        {
            return n * factorial(n - 1);
        }
    }
}
}
4

4 回答 4

3

这些可能是评论而不是答案...

您递归地计算阶乘。迭代计算阶乘可能会更快；与以往一样，您应该对此进行测试，以找出在您的平台上最有效的方法。更糟糕的是，您在阶乘函数的开头测试了一个与 0 相等的 double 值。您使用的函数仅适用于整数，如果您真的想计算实数的阶乘，您应该（最有可能）使用 Gamma Function。

由于从to的定积分等于从 to 的定积分加上从ato的定积分（对于行为合理的函数），您可以将定积分的计算拆分为任意数量的块。baccba < c < b

于 2012-07-25T09:11:08.873 回答
2

阶乘计算效率极低，尤其是当输入值变大时。我还想记住计算 - 为什么一旦有了它就重做呢？

更好的解决方案是使用 gamma 函数来实现它。它也更能抵抗溢出，因为它返回一个双精度值。

我也会担心这个计算的数值错误：

temp = ((Math.Pow(degreeFreedomE, degreeFreedomE / 2) * Math.Pow(degreeFreedomT, degreeFreedomT / 2)) / (factorial(degreeFreedomE / 2 - 1) * factorial(degreeFreedomT / 2 - 1))) * (factorial(((degreeFreedomT + degreeFreedomE) / 2 - 1)))*((Math.Pow(x, degreeFreedomE / 2 - 1)) / (Math.Pow((degreeFreedomT + degreeFreedomE * x), ((degreeFreedomE + degreeFreedomT) / 2))));

您正在乘以和除以可能很大的数字。您希望四舍五入不会杀死您，并且取消会成功

另一件要尝试的事情是利用对数来发挥自己的优势。它们有两个不错的属性：

ln(A*B) = ln(A) + ln(B)

和

ln(A/B) = ln(A) - ln(B)

这将减少那些较大数字的大小，并使计算更不容易出现舍入错误。

于 2012-07-25T09:05:06.383 回答
2

如下更改 for 循环以删除 for 循环内的 if 条件：

for (int i = 1; i < iterations; i++)
{
    x = start + i * dist;
    sumT += integralFunction(x - dist / 2, degreeFreedomT, degreeFreedomE);
    sum += integralFunction(x, degreeFreedomT, degreeFreedomE);
}
x = start + iterations * dist;
sumT += integralFunction(x - dist / 2, degreeFreedomT, degreeFreedomE);

编辑：您也可能想要使用此编译器选项-finline-functions（适用于 gcc 和 icc。检查 C# 中的等效项）。对函数的调用integralFunction （一个简单的 1 行函数）将在编译期间内联，并且可以消除每次迭代的函数调用开销

于 2012-07-25T09:24:13.930 回答
1

您可以使用并行化您的外部 for 循环Parallel.For。请注意，您不能随处使用它。这在很大程度上取决于算法和/或数据。 

public static double integrate(double start, double end, int degreeFreedomT, int degreeFreedomE)
{
    int iterations = 100000;
    double x, dist, sum = 0, sumT = 0;
    dist = (end - start) / iterations;
    Parallel.For(1, iterations, i => {
        x = start + i * dist;
        sumT += integralFunction(x - dist / 2, degreeFreedomT, degreeFreedomE);
        if (i < iterations)
        {
            sum += integralFunction(x, degreeFreedomT, degreeFreedomE);
        }
    });

    sum = (dist / 6) * (integralFunction(start, degreeFreedomT, degreeFreedomE) + integralFunction(end, degreeFreedomT, degreeFreedomE) + 2 * sum + 4 * sumT);
    return sum;
}
于 2012-07-25T09:05:09.620 回答