[Haskell] 分散と標準偏差を求める

Haskellのお勉強

Haskellで平均を求める処理を書くことができたので、そこからさらに分散と標準偏差を求めてみました。

Haskellのコード

main = do
    print $ average numbers           -- 55.0
    print $ variance numbers          -- 600.0
    print $ standartDeviation numbers -- 24.49489742783178
    where numbers = [1..1000]

average :: (Real a) => [a] -> Double
average xs = (realToFrac $ sum xs) / fromIntegral (length xs)

variance :: (Real a) => [a] -> Double
variance xs = average $ foldr (\x acc -> ((realToFrac x - avg)^2):acc) [] xs
  where avg = average xs

standartDeviation :: (Real a) => [a] -> Double
standartDeviation xs = sqrt $ variance xs

コードの解説

上のコードでは平均・分散・標準偏差の3つを求める関数を用意しています。

分散の求め方

分散(variance)を求めるためにまず偏差を求めます。偏差とは個々のデータから平均値を引いた値のことです。

分散は、偏差の2乗をすべて足したものをデータの総数で割ることで求められます。つまり偏差の2乗の平均です。

上の例では引数Real型のリストから、右畳み込みで要素を取り出しながら、偏差の2乗のリストを新たに作成しています。
具体的にはラムダ式の中で、個々の要素をFractional型に変換した上で平均との差を求め、さらに2乗しています。
さらに作成された偏差の2乗リストから平均を求めたものを分散の値としています。

この新しいリストを作成する際に、左畳込みでリストの連結を使うと、データ数が多くなると速度の差が顕著になります。
リストの連結(++)は、先頭への追加(:)に比べてデータ数が多い場合に処理がコストが大きくなるためです。
すごいHaskell本にも、新しいリストを作成する場合は、foldrが良いと書かれていました。

標準偏差の求め方

分散が求められれば標準偏差(standartDeviation)はその正の平方根を求めればよいので、sqrtを使えばおしまいです。

確認

結果が正しく得られているか確認するために.NET向けの数値計算のライブラリで求められる値と比較してみます。使用するのはMath.NET Numericsです。

Nugetからインストールして、MathNet.Numerics.Statisticsをusingに追加してやれば、IEnumerable<double>の拡張メソッドで分散や標準偏差などを求められます。

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

namespace Test
{
    class Program
    {
        static void Main(string[] args)
        {
            var data = new double[] { 20, 30, 40, 50, 60, 70, 80, 90 };
            Console.WriteLine(data.Average());              // 平均 55
            Console.WriteLine(data.Variance());             // 分散 600
            Console.WriteLine(data.StandardDeviation());    // 標準偏差 24.4948974278318
            Console.ReadKey();
        }
    }
}

標準偏差の最終桁が切り上られているようですが、結果が一致していることを確認できました。