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

これらは答えではなくコメントかもしれません...

階乗を再帰的に計算します。階乗を繰り返し計算する方が速い場合があります。いつものように、これをテストして、プラットフォームで最適に機能するものを見つけてください。さらに悪いことに、階乗関数の開始時に double 値が 0 で等しいかどうかをテストします。使用している関数は整数に対してのみ正しいです。実際に実数の階乗を計算したい場合は、(ほとんどの場合) ガンマ関数を使用する必要があります。

からの定積分は からの定積分に to からの定積分を加えたものに等しいのでa(適切にb動作する関数の場合)、定積分の計算を必要な数のチャンクに分割できます。accba < c < b

于 2012-07-25T09:11:08.873 に答える
2

階乗計算は、特に入力値が大きくなると、非常に効率が悪くなります。私も計算をメモしたいと思います。

より良い解決策は、ガンマ関数を使用して実装することです。double 値を返すため、オーバーフローに対する耐性も高くなります。

この計算では、数値エラーも心配されます。

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))));

大きい可能性のある数値を乗算およ​​び除算しています。ラウンドオフがあなたを殺さず、キャンセルがうまくいくことを望んでいます

試してみるべきもう 1 つのことは、対数を有利に利用することです。これらには 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 に答える