1. 幾何平均 gmean
gmean(a, axis=0, dtype=None, weights=None)
from scipy.stats import gmean
x = [1.2, 3.1, 4.6, 2.2, 5.7]
gmean(x)
2.9263054868319522
gmean 関数を使わずに計算すると以下のようになる。
import numpy as np
np.exp(np.mean(np.log(x)))
2.9263054868319522
np.prod(x)**(1/len(x)) # 非推奨
2.9263054868319527
以下のようなテストデータ（平均値=50，標準偏差=10)で幾何平均を求めてみる。
from scipy.stats import norm
z = norm.rvs(loc=50, scale=10, size=200)
import matplotlib.pyplot as plt
plt.hist(z, bins=30);
サンプルサイズ 200 で，教科書的な定義式通りの np.prod(z)**(1/len(z)) は途中でオーバーフローしてしまい inf を返す。
gmean(z), np.exp(np.mean(np.log(z))), np.prod(z)**(1/len(z)), np.prod(z)
(49.19589704885539, 49.19589704885539, inf, inf)
2. 調和平均 hmean
hmean(a, axis=0, dtype=None)
from scipy.stats import hmean
x = [1.2, 3.1, 4.6, 2.2, 5.7]
hmean(x)
2.4958950838709315
hmean 関数を使わずに計算すると以下のようになる。
1/np.mean(1/np.array(x))
2.4958950838709315
`np.array` を使ったのは，`x` がリストのとき `1/x` はエラーになるので，numpy.ndarray に変換するためである。
どうせ `numpy` を使うなら，逆数を返す `np.reciprocal` を使えば，統一性がある。
np.reciprocal(np.mean(np.reciprocal(x)))
2.4958950838709315