文系プログラマが数式アレルギーを治す為に、中学〜高校レベルの知識でもわかる数式を、プログラムのコードに直す作業を黙々と続ける企画。
ソースコードはJuliaを利用。2012年生まれ、もうすぐ3歳という幼さにして完成度たけーなこのプログラミング言語、Juliaを利用。
前回、愚直に乗算で幾何平均を出したけど、こんなんちょっと大きい桁扱わせたら簡単に溢れるじゃないかということで、対数を使って加算で計算する。
式はWikipediaの幾何平均を参考にした。
http://ja.wikipedia.org/wiki/%E5%B9%BE%E4%BD%95%E5%B9%B3%E5%9D%87
こんな感じ。
# 自分用メモ
$ \exp( \frac {1} {n} \sum_{i=1}^n \ln {x_i} ) $
Σ（大文字のシグマ）は集合の和を表す。例えば「2, 5, 8, 3, 6」の集合だったら、「2+5+8+3+6」になる。
lnは自然対数。eを底とする対数。プログラムだとlog(x)と書くとeが対数とされる場合が多い。
expはeを底とする指数関数。logと逆の意味になるので、log(x)をexpすると元の値に戻し、exp(x)をlogしても元に戻る。
exp( log( pi ) ) #=> 3.141592653589793 log( exp( pi ) ) #=> 3.141592653589793
算術平均で出てきた数式と見比べると、本式はlogの算術平均のexpということがわかる。対数を使うと乗算を加算できるのでこれで幾何平均が取れる。
ということで、これをJuliaのコードに直すと、下記のようになる。
x = [3.0, 2.5, 1.4, 6.2, 3.3, 2.1, 5.0] x_sum = 0.0 for i = x x_sum += log(i) end n = length( x ) exp( (1/n) * x_sum ) #=> 3.013287735358114
これで前回やったscipyの値（3.013287735358115）と一致……してないね。
なんでや。scipyのソースを要約すると、式はこうか。
using PyCall @pyimport numpy as np np.exp( mean( np.log(x) ) ) #=> 3.013287735358115
これをできるだけ近い感じでjuliaに直すと、こんな感じだろうか。
exp( mean( map( a -> log(a), x ) ) ) #=> 3.013287735358114
これは浮動小数点の扱いに差があるのか。この先を調べるのは時間がかかりそうなので、ここで撤退する。