Is my integration algorithm correct? Also can we improve this on parallel processing and efficiency? The code is in C#

Question

How do i better calculate the definite integral? I am using a function to integrate and another to find the factorial recursively.

I'l like to better the algorithm or the efficiency or even the accuracy for that matter.

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

Answer 1

change the for loop as below to remove the if condition inside the for loop:

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

EDIT: Also you might want to use this compiler option -finline-functions(works in gcc and icc. check for equivalent in C#). The call to the function integralFunction (a simple 1 line function) will be inlined during compilation and function call overhead for each iteration can be removed

Answer 2

These might be comments rather than an answer ...

You compute the factorial recursively. It might be faster to compute the factorial iteratively; as ever you should test this to find out what works best on your platform. Worse, you test a double value for equality with 0 at the start of your factorial function. The function you are using is correct only for integers, if you really want to compute factorials for real numbers you should be using (most likely) the Gamma Function.

Since the definite integral from a to b is equal to the definite integral from a to c plus the definite integral from c to b (for a < c < b and for reasonably well-behaved functions) you can split the computation of the definite integral into as many chunks as you wish.

Answer 3

The factorial calculation is terribly inefficient, especially as the input value gets larger. I'd also want to memoize the calculation - why redo it once you have it?

A better solution is to implement it using the gamma function. It'll be more resistant to overflow as well, because it returns a double value.

I'd also be worried about numerical errors with this calculation:

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

You're multiplying and dividing numbers that may be large. You're hoping that roundoff doesn't kill you and that cancellations work out

Another thing to try is to exploit logarithms to your advantage. They have two nice properties:

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

and

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

This will keep the size of those larger numbers down and make the calculation less prone to round off errors.

Answer 4

You can parallelise your outer for loop using Parallel.For. Please note that you cannot just use this anywhere and everywhere. It depends a lot on the algorithm and/or data.

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