普段NumPyやPandas、SQLなどを使っているとそもそも関数が用意されていて意識しなくとも使えますが、久々の復習も兼ねて自分で書いたりして理解を深めます。言語はPythonを使います（Pythonビルトインのパッケージもぼちぼち使っていきます）。
平均は触れることが少ないのでスキップして、中央値、最頻値、分散、標準偏差、相関係数について計算します。
※著者は理系出身などではなくデザインの学校卒の人間ですので知識や理解が粗い点などがご容赦ください。
使う環境
- Python 3.9.0 (+ Jupyter)
- NumPy 1.20.3 （検算などでのみ使用します）
中央値
中央値（median）はデータ群を値の大きさでソートした状態で中央に位置している値です。シンプルな計算ですが、データ群の件数が奇数偶数で計算が変わる点には注意が必要です。
奇数の場合にはそのまま中央の位置に位置する値を取ればOKですが、偶数の場合は中央に位置する値が2つとなるためその2つを合算して2で割る必要があります。
データ群の件数をNとすると、奇数の場合には中央の値の位置は以下のようになります。
\frac{(N + 1)}{2}
一方で偶数の場合には以下の2つの位置を取得する必要があります。
\tag{1}
\frac{N}{2}
\tag{2}
\frac{N}{2} + 1
※Python上でリストなどで扱う場合にはインデックスが0からスタートするのでこれらの値からさらに1を減算するため、$\frac{N}{2} - 1$と$\frac{N}{2}$の2つになります。
from typing import Union, List
def calculate_median(numbers: List[Union[int, float]]) -> Union[int, float]:
"""
指定されたデータ群のリストの中央値を取得します。
Parameters
----------
numbers : list
対象の数値のリスト。
Returns
-------
median : int or float
計算結果の中央値。
"""
# データの件数を取得しています。
N: int = len(numbers)
# 中央値の計算では事前にソートが必要になります。
numbers = sorted(numbers)
# 奇数・偶数判定を行っています。
if N % 2 != 0:
# 奇数の場合の計算です。シンプルに1つの位置の値を取れば計算できます。
index: int = int((N + 1) / 2) - 1
median: Union[int, float] = numbers[index]
return median
# 偶数の場合の計算です。2つの値の参照が必要になります。
index_1: int = int(N / 2) - 1
index_2: int = int(N / 2)
median = (numbers[index_1] + numbers[index_2]) / 2
return median
検算
少し動かしたりNumPy側での計算結果との比較などを行ってみます。
まずは奇数の件数を持つリストを指定してみます。
numbers: List[Union[int, float]] = [31, 21, 39, 10, 23, 41, 53, 18, 9, 18, 25]
median: Union[int, float] = calculate_median(numbers=numbers)
print(median)
23
続いて偶数の件数を持つリストを指定しています。こちらは2で除算されるので結果がfloatになっています。
numbers = [31, 21, 39, 10, 23, 41, 53, 18, 9, 18, 25, 12]
median = calculate_median(numbers=numbers)
print(median)
22.0
NumPyでも検算で計算してみます。
import numpy as np
median = np.median([31, 21, 39, 10, 23, 41, 53, 18, 9, 18, 25])
print('奇数の場合:', median)
median = np.median([31, 21, 39, 10, 23, 41, 53, 18, 9, 18, 25, 12])
print('偶数の場合:', median)
奇数の場合: 23.0
偶数の場合: 22.0
一致しました。大丈夫そうです。
最頻値
今度は最頻値（mode）を計算してみます。データ群の中で一番頻度が多く出てくる値が対象になります。
各値が何回出てくるのか・・・という計算をしたい場合にはビルトインを使う場合にはcollectionsパッケージのCounterクラスが便利です。
コンストラクタにリストなどのiterableオブジェクトを渡すと各値ごとの出現数をカウントしてくれます。
most_commonメソッドを使うと2つの値を格納したタプルのリストが返ってきます。タプルの先頭の値は対象の値、2つ目の値は頻度（出現数）となります。頻度が多い値から順番にリストに格納されます。
from collections import Counter
counter: Counter = Counter([1, 3, 1, 1, 2, 3])
print(counter.most_common())
[(1, 3), (3, 2), (2, 1)]
上位N件のみの値のみ欲しければNの値をmost_commonメソッドの引数に指定します。今回は最頻値で先頭の値だけあれば良いので1を指定します。
print(counter.most_common(1))
[(1, 3)]
最頻値では出現回数ではなく対象の値が欲しいので返却値に対して[0][0]というインデックスの指定をします（結果は一番リスト内で多い1となります）。
print(counter.most_common(1)[0][0])
1
関数にすると以下のような感じになります。
from typing import Iterable, TypeVar
from collections import Counter
T = TypeVar('T')
def calculate_mode(iterable: Iterable[T]) -> T:
"""
指定されたデータ群のリスト内の最頻値を取得します。
Parameters
----------
iterable : Iterable
対象のリストなどのデータ群。
Returns
-------
mode : Any
最頻値の値。
"""
counter: Counter = Counter(iterable)
mode: T = counter.most_common(1)[0][0]
return mode
print(calculate_mode([1, 3, 1, 1, 2, 3]))
1
分散
続いて分散（variance）を計算していきます。分散は以下のような計算になります。
分散とは、データの散らばりの度合いを表す値です。分散を求めるには、偏差（それぞれの数値と平均値の差）を二乗し、平均を取ります。
分散の意味と求め方、分散公式の使い方
分散（variance）を求める計算を式で表すと以下のようになります。分母のnはデータの件数、$x_i$はデータ群のi番目の値、$x_{\textrm{mean}}$はデータ群の平均値です。
\textrm{variance} =
\frac{\sum_{n=1}^n (x_i - x_{\textrm{mean}})^2}
{n}
Pythonで書いていきましょう。
from typing import List, Union
def calculate_variance(numbers: List[Union[int, float]]) -> float:
"""
指定されたデータ群の分散を取得します。
Parameters
----------
numbers : list
対象のデータ群のリスト。
Returns
-------
variance : float
計算後の分散。
"""
x_mean: float = sum(numbers) / len(numbers)
squared_diff_sum: float = 0
for x_i in numbers:
diff: float = x_i - x_mean
squared_diff_sum += diff ** 2
n: int = len(numbers)
variance: float = squared_diff_sum / n
return variance
まずはx_mean: float = sum(numbers) / len(numbers)という部分で平均を出しています。式の$x_{\textrm{mean}}$という部分が該当します。
diff: float = x_i - x_meanという部分は式の偏差の$x_i - x_{\textrm{mean}}$部分に該当します。
2乗計算はdiff ** 2と書けるのと、式の$\sum_{n=1}^n$の合計を出す計算のためにsquared_diff_sum += diff ** 2と書いて合計を計算しています。
分母のnはlen(numbers)でリストの長さで取っています。
検算
実際に実行してみます。
calculate_variance(numbers=[100, 120, 90, 250, 160, 180, 200, 80, 120, 170])
2661.0
NumPyでも分散計算をしてみて比較してみます。NumPyではvar関数で分散を取れます。
np.var([100, 120, 90, 250, 160, 180, 200, 80, 120, 170])
2661.0
NumPy側と値が一致しました。実装分が合っていそうです。
標準偏差
標準偏差（standard deviation）は分散の値の平方根を取ることで計算ができます。式にすると以下のようになります。
\textrm{standard deviation} = \sqrt{\textrm{variance}}\
分散計算の関数はもう既に前節で触れているので計算は非常にシンプルになります。
Pythonビルトインでの平方根の計算はmathモジュールのsqrt関数を使うかもしくは** 0.5と計算します。どちらも同じ結果となります。
import math
print(math.sqrt(2661.0))
51.58488150611572
print(2661.0 ** 0.5)
51.58488150611572
関数にすると以下のようになります。分散計算の関数部分を除けば非常にシンプルです。
from typing import List, Union
def calculate_variance(numbers: List[Union[int, float]]) -> float:
"""
指定されたデータ群の分散を取得します。
Parameters
----------
numbers : list
対象のデータ群のリスト。
Returns
-------
variance : float
計算後の分散。
"""
x_mean: float = sum(numbers) / len(numbers)
squared_diff_sum: float = 0
for x_i in numbers:
diff: float = x_i - x_mean
squared_diff_sum += diff ** 2
n: int = len(numbers)
variance: float = squared_diff_sum / n
return variance
def calculate_standard_deviation(numbers: List[Union[int, float]]) -> float:
"""
指定されたデータ群の標準偏差を取得します。
Parameters
----------
numbers : list
対象のデータ群のリスト。
Returns
-------
standard_deviation : float
計算後の標準偏差。
"""
variance: float = calculate_variance(numbers=numbers)
standard_deviation: float = variance ** 0.5
return standard_deviation
検算
追加した関数で計算してみます。
calculate_standard_deviation(
numbers=[100, 120, 90, 250, 160, 180, 200, 80, 120, 170])
51.58488150611572
NumPy側でも計算して比較してみます。NumPy側では標準偏差はstd関数で取れます。
import numpy as np
np.std([100, 120, 90, 250, 160, 180, 200, 80, 120, 170])
51.58488150611572
NumPy側の値と一致しています。大丈夫そうです。
相関係数
最後に相関係数（correlation coefficient）の計算を書いていきます（本記事ではピアソンの積率相関係数（Pearson correlation coefficint）で進めます）。
1つ目のデータ群を$x$, 2つ目のデータ群を$y$、各データの長さを$n$とすると相関係数の計算は以下のようになります。
\textrm{correlation coefficient} =
\frac
{n\sum xy - \sum x\sum y}
{\sqrt{(n\sum x^2 - (\sum x)^2)(n\sum y^2 - (\sum y)^2)}}
若干式は長めですが、分母部分は同じような計算をxとyそれぞれでやっている都合で長くなっているのと、やっていることとしては総和や2乗、平方根や掛け算や引き算・・・といった具合なので内容としては簡単な計算となります。
分母の部分の$\sum x^2$と$(\sum x)^2$部分は括弧の有る無しのみで記述が似ていますが計算内容が異なりますので注意します。前者は個別の各値を2乗してからの総和となり、後者は総和してからの2乗となります。y側も同様です。
Pythonで書いていきます。若干無駄が発生しますが今回は分かりやすさのために分子と分母の計算をそれぞれ関数を分けています。
まずは分子（numerator）部分の計算から対応します。
from typing import Union, List
def calculate_correlation_coefficient_numerator(
x: List[Union[int, float]],
y: List[Union[int, float]]) -> float:
"""
指定された2つのデータ群の相関係数計算の分子側の値を取得します。
Parameters
----------
x : list
1つ目の対象のデータ群のリスト。
y : list
2つ目の対象のデータ群のリスト。
Returns
-------
numerator : float
相関係数計算の分子側の計算結果。
"""
n: int = len(x)
xy_prods: List[float] = []
for x_i, y_i in zip(x, y):
xy_prods.append(x_i * y_i)
numerator: float = n * sum(xy_prods) - sum(x) * sum(y)
return numerator
必要な計算式部分は$n\sum xy - \sum x\sum y$です。
nはデータ群の件数を取れば良いのでlen(x)で取っています（基本的にxとyが同じ件数となる想定なのでこれはxでもyでも変わりません）。
$\sum xy$部分の計算用に個別のxとy同士を掛けた値が必要になるのでループで回してxy_prodsというリストにx_i * y_iという計算値を格納しています。最終的にn * sum(xy_prods) - sum(x) * sum(y)という数式通りの値を計算して分子側は完了です。
続いて分母（denominator）側の計算を対応してきます。
def calculate_correlation_coefficient_denominator(
x: List[Union[int, float]], y: List[Union[int, float]]) -> float:
"""
指定された2つのデータ群の相関係数計算の分母側の値を取得します。
Parameters
----------
x : list
1つ目の対象のデータ群のリスト。
y : list
2つ目の対象のデータ群のリスト。
Returns
-------
denominator : float
相関係数計算の分母側の計算結果。
"""
n: int = len(x)
squared_x_list: List[Union[int, float]] = []
for xi in x:
squared_x_list.append(xi ** 2)
squared_x_sum: Union[int, float] = sum(squared_x_list)
x_term: Union[int, float] = n * squared_x_sum - sum(x) ** 2
squared_y_list: List[Union[int, float]] = []
for yi in y:
squared_y_list.append(yi ** 2)
squared_y_sum: Union[int, float] = sum(squared_y_list)
y_term: Union[int, float] = n * squared_y_sum - sum(y) ** 2
denominator: Union[int, float] = (x_term * y_term) ** 0.5
return denominator
必要な式は$\sqrt{(n\sum x^2 - (\sum x)^2)(n\sum y^2 - (\sum y)^2)}$です。
$n\sum x^2 - (\sum x)^2$部分の計算をx_term、$n\sum y^2 - (\sum y)^2$部分をy_termという変数名にしています。
括弧無しの$\sum x^2$といった部分では個別の値を先に2乗してから総和を取る必要があるため、ループで回してsquared_x_list.append(xi ** 2)という記述でリスト内に個別の値の2乗値を格納しています。
括弧付きの$(\sum x)^2$の方の計算はループなどは回さずにxのリストの総和を2乗すれば良い形となります（sum(x) ** 2という部分が該当します）。
y側も同じように計算したら、最後にx側の計算（x_term）とy側の計算（y_term）を掛けた後に平方を取って分母側は完成です。x_termなどを反映して式で表すと$\sqrt{(\textrm{x_term})(\textrm{y_term})}$といったようになります。平方は** 0.5という計算で取れるのでコード上では(x_term * y_term) ** 0.5という部分が該当します。掛け算などとは異なり平方を反映する部分には括弧が無いと計算結果が変わってしまうので注意します。
最後に分子と分母による割り算を行うための関数を追加して完成です。
def calculate_correlation_coefficient(
x: List[Union[int, float]], y: List[Union[int, float]]) -> float:
"""
指定された2つのデータ群の相関係数を取得します。
Parameters
----------
x : list
1つ目の対象のデータ群のリスト。
y : list
2つ目の対象のデータ群のリスト。
Returns
-------
correlation_coefficient : float
計算された相関係数。
"""
numerator: float = calculate_correlation_coefficient_numerator(x=x, y=y)
denominator: float = calculate_correlation_coefficient_denominator(
x=x, y=y)
correlation_coefficient: float = numerator / denominator
return correlation_coefficient
検算
実行してみて相関係数を取得してみます。
calculate_correlation_coefficient(
x=[10, 8, 9, 6, 3, 5, 2], y=[22, 19, 18, 15, 16, 14, 13])
0.8631908291055478
約0.86ということで強い正の相関値が出ています。
NumPy側でも計算してみます。NumPyではcorrcoef関数を使います。3つ以上のデータ群を指定することがあるため結果は行列となりますが、今回は2つの値でしか使わないので一部の値のみ参照します（1ではない部分の値を参照します）。
import numpy as np
np.corrcoef([10, 8, 9, 6, 3, 5, 2], [22, 19, 18, 15, 16, 14, 13])
array([[1. , 0.86319083],
[0.86319083, 1. ]])
NumPy側も約0.86と出ており一致しています。大丈夫そうです。