Lispで統計の話(平均・中央値・分散) - たくさん寝太郎の寝床

こんにちは、たくさん寝太郎です。

先日の記事でランダムウォークについて考えたのですが、ついでに統計量についても気になったので考えてみました。

以下の内容はこの本を参考にして書いています。

新統計入門

作者: 小寺平治
出版社/メーカー: 裳華房
発売日: 1996/10/01
メディア: 単行本
この商品を含むブログを見る

相加平均・中央値

数のリスト $(x_1,x_2,...,x_N)$ が与えられた時、相加平均は次のように定義されます。

$\bar{x} = {\displaystyle \frac{1}{N}\sum_{k=1}^{N} x_k}$

まずリストの総和を求める関数を定義しましょう。

(defun sum (num_list)
  (if num_list
      (+ (car num_list)
	 (sum (cdr num_list)))
      0))

以前も話した通り、Lispでは空のリストを偽として扱います。リストの先頭+残りのリストの総和という再帰処理にしてあげれば総和が定義できます。

(追記 2018/09/14)
applyコマンドを用いれば再帰的な処理をしなくても総和が計算できます。

(apply #'+ num_list)

これにより相加平均が定義できます。

(defun ave (num_list)
  (/ (sum num_list)
     (length num_list)))

中央値はN個のデータをソートした時に中央にくる値です。*1

(defun med (num_list)
  (let* ((sort_list (sort num_list #'<))
        (N (length num_list))
        (m (if (oddp N)
               (/ (1- N) 2)
               (/ N 2))))
        (if (oddp N)
            (nth m sort_list)
            (/ (+ (nth (1- m) sort_list)
		  (nth m sort_list))
	       2)
            )))

こんな感じでしょうか。もっと楽な書き方もありそうですが...。
注意点ですが、local変数定義でlet*を使っているのはNを定義した後mの定義でNを用いているからです。letを使うと「IF: variable N has no value」とerrorが出てきます。

乱数(0~99)を1000000個生成してリストを作成し、中央値を求めるのにどれくらい時間がかかるか試してみました。

$time clisp stat.lisp
50.0
clisp stat.lisp  10.59s user 4.30s system 98% cpu 15.041 total

ちょっと遅いですね...

(余談)
Pythonで同じく乱数を1000000個生成して中央値を求めると次のようになりました。

$time python stat.py
50.0
python stat.py  2.29s user 0.07s system 95% cpu 2.468 total

分散

数のリスト $(x_1,x_2,...,x_N)$ が与えられた時、平均を $\bar{x}$ とすると分散 $\sigma^2(x)$ は次のように定義されます。

$\sigma^2(x) = {\displaystyle \frac{1}{N}\sum_{k=1}^N(x_k-\bar{x})^2}$

この式は次と同値です。

$\sigma^2(x) = {\displaystyle \frac{1}{N}\sum_{k=1}^{N} x_k^2} - {\displaystyle (\frac{1}{N}\sum_{k=1}^{N} x_k)^2}$

(defun var (num_list)
  (- (ave (map 'list #'* num_list num_list))
     (expt (ave num_list) 2)))

また、分散の平方根を標準偏差と呼びます。

(defun sd (sigma2)
  (expt sigma2 1/2))

今回はこの辺りで終わりです。
次回は相関係数あたりの実装をしようと思います。

*1:厳密な定義は以下参照:中央値 - Wikipedia