.NETでライブラリを使わずLINQで標準偏差(StdDev)を計算する方法

この記事は公開から2年以上経過しています。

.NETには標準偏差を計算する標準ライブラリ関数がないので、C#/VB.NETのLINQで標準偏差(StdDev/Standard Deviation/STDEV.P)を計算する簡単な拡張メソッドを自作してみたので紹介します。

Math.NET Numerics等、社外ライブラリを使えない場合に役立つかも知れません。

サンプルソースコード

実装方法によって浮動小数点による小数点誤差が生じるため、有効桁数に注意が必要です。

C#

static class Extensions
{
    public static double StdDev(this IEnumerable<double> src)
    {
        IList<double> lst = src is IList<double> ? (IList<double>)src : src.ToArray();
        int n = lst.Count();
        double ave = lst.Average();
        double s = 0;
        for (int i = 0; i < n; ++i)
            s = s + (lst[i] - ave) * (lst[i] - ave);
        return Math.Sqrt(s / n);
    }
}

VB.NET

Module Extensions

    <Extension()>
    Public Function StdDev(src As IEnumerable(Of Double)) As Double
        Dim lst As IList(Of Double) = If(TypeOf src Is IList(Of Double), CType(src, IList(Of Double)), src.ToArray())
        Dim n As Integer = lst.Count()
        Dim ave As Double = lst.Average()
        Dim s As Double = 0
        For i As Integer = 0 To n - 1
            s = s + (lst(i) - ave) * (lst(i) - ave)
        Next
        Return Math.Sqrt(s / n)
    End Function

End Module

GitHub Gist

実装は標準偏差の計算式

の通りです。

標本に対する標準偏差(STDEV.S)の場合はメソッド末尾をMath.Sqrt(s / (n - 1))に変更してください。

標本数が1以下の場合は結果がNaNになります。

検証

以下のようなコードで10個の乱数の実数に対してMath.NET NumericsのPopulationStandardDeviation()と本メソッドを使って標準偏差を求めてみた結果は以下の通りです。

テストコード(C#)

using System;
using System.Collections.Generic;
using System.Linq;
using MathNet.Numerics.Statistics;

private static void Main(string[] args)
{
    var rnd = new Random();

    for (var i = 0; i < 10; ++i)
    {
        var src = Enumerable.Range(0, 10).Select(_ => rnd.Next(-10000000, 10000001) / 10000d).ToArray();

        var stddev1 = src.StdDev();
        var stddev2 = Statistics.PopulationStandardDeviation(src);
        Console.WriteLine($"[{string.Join(",", src)}]");
        Console.WriteLine($"stddev1:{stddev1}");
        Console.WriteLine($"stddev2:{stddev2}");
    }
}

出力結果：

stddev1:378.961251135221
stddev2:378.961251135221
[-501.9189,335.5237,400.4169,-719.991,-447.1985,-435.2201,-121.6716,56.6706,124.0449,-49.3778]
stddev1:358.308644842964
stddev2:358.308644842964
[-785.5789,126.4972,-626.6209,541.5176,9.3645,408.6149,-72.7582,-652.9837,-823.9302,910.5459]
stddev1:577.28107311123
stddev2:577.28107311123
[-191.3754,533.7034,220.7657,604.4507,-145.7344,-370.0983,-276.7076,-613.4914,41.6237,302.2906]
stddev1:379.303667358862
stddev2:379.303667358862
[-257.0222,45.2968,48.0518,913.5202,692.7316,324.156,952.9933,-870.2568,737.1815,792.1251]
stddev1:564.333481740432
stddev2:564.333481740432
[-700.633,267.3772,886.9665,-550.5951,270.802,338.9406,-960.9497,467.3376,-302.2009,880.4755]
stddev1:617.506624032449
stddev2:617.506624032449
[719.3618,-788.5392,-164.0506,29.0006,-899.0651,-35.6223,-885.242,795.9661,425.7433,922.8998]
stddev1:665.235451315808
stddev2:665.235451315808
[841.4875,-508.1802,-984.9531,821.7542,-792.3556,169.1242,-921.5547,-540.1411,848.7859,685.2084]
stddev1:747.305298597112
stddev2:747.305298597112
[455.6216,-202.5014,-390.1652,-616.0554,611.1788,541.2474,646.7368,-228.6392,-666.3513,515.325]
stddev1:508.188847331742
stddev2:508.188847331742
[-409.1487,255.9547,688.7094,161.1482,582.9287,-468.483,-790.6526,-980.9489,812.5853,-271.6461]
stddev1:599.696866548962
stddev2:599.696866548962

ちなみに私のLinuxマシン(Xeon E3-1270 v6)上のQEMU KVM仮想環境上のWindows 10でデバッガ経由で10万個の数値を単純ループで連続1000回計算するのに要した時間は

名称	計算時間(ms)
Math.NET Numerics PopulationStandardDeviation()	520
本拡張メソッド	1,480

と、本拡張メソッドの計算速度は、Math.NETの約1/3でした。
※内部処理でMath.Pow()を使うと更に遅くなり約1/10。
パフォーマンスを求めるのであればMath.NETを使うか、ロジックを転用したほうが良いかと思います。
(Welfordアルゴリズムを使えば1回のループで処理可能です。)

参考ウェブサイトなど

Math.NET
Math.NET Numerics

以上です。