Parzen Window Density Estimation in C# Is my implementation correct (it is matching this result)? How can I improve this code?

public class Parzen
{
    public static double Variance(IEnumerable&lt;double&gt; source)
    {
        double avg = source.Average();

        double d = source.Aggregate(0.0, (total, next) =&gt; total += Math.Pow(next - avg, 2));

        return d / (source.Count() - 1);
    }

    private static double GausianFunction(double x1, double x, double sigma2)
    {
        double up = Math.Pow(x1 - x, 2);

        double dw = 2 * sigma2;

        double firstPart = 1 / Math.Sqrt(2 * Math.PI * sigma2);

        double ret = firstPart * Math.Exp((-1) * (up / dw));

        return ret;
    }

    public static double PrazenWindowDensity(double [][] Xn, double x, double sigma2)
    {
        double gauss = 0;

        foreach(double [] arr in Xn)
        {
            foreach (double item in arr)
            {
                double xx = GausianFunction(item, x, sigma2);
                gauss += xx;
            }
        }

        return gauss / Xn.Sum(xy =&gt; xy.Length); ;
    }
}

public class Program
{
    static void Main(string[] args)
    {
        double[][] mat = 
        {
            new double[]{2, 2.5}, 
            new double[]{3, 1},
            new double[]{6}
        };

        //double sigma2 = Parzen.Variance(mat.Cast&lt;double&gt;().ToList());

        double den = Parzen.PrazenWindowDensity(mat, 3, 1);

        string str = string.Empty;
    }
}
I assume your implementation is correct then I go directly to (2).

Even assuming you do not worry too much about speed you should replace foreach with for in PrazenWindowDensity(), relative speed-up is consistent:

for (int m=0; m &lt; xn.Length; ++m)
{
    for (int n=0; n &lt; xn[m].Length; ++n)
    {
        gauss += GausianFunction(xn[m][n], x, sigma2);
    }
}


Few notes: I renamed Xn to xn because parameters should be camelCase. I dropped local variable xx because it's not helpful here. 

In the same function you're using LINQ to sum lengths of each array, it's deadly slow and it may be replaced with a simple ad-hoc function:

static int SumLengths(double[][] data)
{
    int sum = 0;
    for (int i=0; i &lt; data.Length; ++i)
        sum += data[i].Length;
}


To have this:

return gauss / SumLengths(xn);




Optimization of GausianFunction() might be more controversial. I observed that often (but not always) C# compiler + JIT compiler will emit better assembly code when an expression is computed without intermediate results (if possible all values stay in 80 bit FPU registers, better accuracy and no need to convert back/to 64 bit double values). I can't say if this is your case but you may give it a chance (let me also use using static Math; I think it helps clarity when I have to work with math formulas):

private static double GausianFunction(double x1, double x, double sigma2)
    =&gt; (1 / Sqrt(2 * PI * sigma2)) * Exp(Pow(x1 - x, 2) / (-2 * sigma2));




Variance() function is actually unused then I'll go quickly through it but: if possible use double[] instead of IEnumerable&lt;double&gt; and avoid LINQ when you can write an ad-hoc function, relative speed-up is huge in this case (of course unless actual performance are enough for your task.)