幾何平均

概要

文系プログラマが数式アレルギーを治す為に、中学〜高校レベルの知識でもわかる数式を、プログラムのコードに直す作業を黙々と続ける企画。

ソースコードはJuliaを利用。Rほどではないけど名前が検索しづらいJuliaを利用。

今回のテーマは幾何平均。相乗平均。geometric mean。scipyだとgmeanという関数名で用意されてる。

数式

式はWikipediaの幾何平均を参考にした。

http://ja.wikipedia.org/wiki/%E5%B9%BE%E4%BD%95%E5%B9%B3%E5%9D%87

こんな感じ。

${\left({\displaystyle \prod _{i=1}^n}{x}_i\right)}^{\genfrac{}{}{0.1ex}{}{1}{n}}$

なんかうちのブラウザだと 1/n のところがちょっと小さかったのでフォントサイズを大きめにした。

# 自分用メモ
$ ( \prod_{i=1}^{n} x_i ) ^ \frac {1} {n} $

数式の解説

Π は 大文字のパイ。集合の積を表す。例えば「2, 5, 8, 3, 6」の集合だったら、「2×5×8×3×6」になる。

Πが含まれるカッコの中だけをコードに直すと、こんな感じ。

x = [3.0, 2.5, 1.4, 6.2, 3.3, 2.1, 5.0]
x_prod = 1
for i = x
  x_prod *= i
end

x_prod
  #=> 132.68911764705882

これを要素数分の1乗したものが、相乗平均になるということらしい。

コードに直す

ということで、これをJuliaのコードに直すと、下記のようになる。

x = [3.0, 2.5, 1.4, 6.2, 3.3, 2.1, 5.0]
x_prod = 1.0
for i = x
  x_prod *= i
end

n = length( x )
x_prod ^ (1/n)
  #=> 3.013287735358114

3.0132くらいの値になった。前回、同じ値を相加平均した時は、3.357くらいの値になっていた。

この計算結果が正しいか、scipyのgmean（幾何平均取ってくれる関数）で確認しておく。

Juliaでは幾何平均を取る関数が見当たらなかったので、PyCallを使ってSciPyを呼び出す。

x = [3.0, 2.5, 1.4, 6.2, 3.3, 2.1, 5.0] 

Pkg.add( "PyCall" )
using PyCall
@pyimport scipy.stats.mstats as mstats
mstats.gmean( x )
  #=> 3.013287735358115

若干誤差が出ているけど(最後の桁が1違う)、概ね正しい数字になっているようだ。

このズレはgmeanと計算方法が違う（対数を利用して計算している）せいだろうか。

今回の計算方法だと、要素数が多かったり桁が大きいと簡単に64bitの精度を超えてしまうので、次回はscipyを見習って対数を使って計算してみることにする。