分散と共分散をPythonで計算する

分散、標準偏差の計算

Pythonを使い分散や共分散を計算します。分散や共分散についてはNumPyモジュールの関数を使うと簡単に求めることができますが、それだけでは身もふたもないので、はじめに定義に従って計算してみます。まずは使う変数をまとめておきます

用語	変数	意　　　　　味
$x$のリスト	list_x	分散などを計算する確率変数　1, 4, 3, 5, 2をリストで定義
$x$のndarray による配列	arr_x	list_xと同じデータをNumPyで計算する場合に使用
平均	mean_x	平均はaverageよりmeanを使うことが多い
偏差平方和	dev_x	偏差平方和はdeviation sum of squares
$x$の2乗の合計	square_x	$x$の2乗を合計するための変数
分散	var_x	分散はvariance
標準偏差	std_x	標準偏差はstandard deviation

定義通りの分散の計算

for文による偏差平方和、分散、標準偏差の計算

list_xを定義して、分散を計算します。最も基本的な方法はfor 文を使った計算です。リストの要素一つ一つについて、for 文を使って偏差の2乗を足し上げていきます。このことにより偏差平方和を求めることができるので、データ数で割って分散、さらに平方根を取ることにより標準偏差を計算します。定義による方法は、各データから平均を一度計算し、さらに各データとその平均の差を求める必要があるので煩わしくなります。

for文を使い分散を計算

list_x = list([1, 4, 3, 5, 2])
mean_x = sum(list_x)/len(list_x)
dev_x = 0
for i in list_x:
dev_x += (i-mean_x)**2
var_x = dev_x/len(list_x)
std_x = var_x**0.5
print("xの偏差平方和:{:.1f} 分散:{:.1f} 標準偏差:{:.1f}".format(dev_x, var_x, std_x))

結果は次の通りです。

xの偏差平方和:10.0  分散:2.0  標準偏差:1.4

リスト内包表記による方法

リスト内包表記を使うと、同じ考え方で簡潔に表記することができます。なお今後は、偏差平方和を求める式のみとします。分散、標準偏差は#1と同じ方法で計算することができます。

リスト内包表記を使い偏差平方和を計算

mean_x = sum(list_x)/len(list_x)
sum([(i-mean_x)**2 for i in list_x])

ラムダ式とmap関数による方法

ラムダ式による方法です。ラムダ式はmap関数と組み合わせてmap(ラムダ式,リスト)とすることで、リストの要素一つ一つにつきラムダ式で定義した計算結果を返します。ただし、map関数の結果はmapオブジェクトで返され、実体としてのリストではないので、通常はlist関数を使いリストに変換します。

ラムダ式～map関数を使い偏差平方和を計算

mean_x = sum(list_x)/len(list_x)
sum(list(map(lambda i: (i-mean_x)**2, list_x)))

なお、今回はmap関数の結果をsum関数で合計しているので、次のようにlistを省略しても問題なく計算することができます。

sum(map(lambda i: (i-mean_x)**2, list_x))

NumPyモジュールによる方法

NumPyモジュールではブロードキャストという機能があり、for文、リスト内包表記、ラムダ式を使わなくても、sum((arr_x-mean_x)**2)だけで、リストの要素一つ一つについて偏差の2乗を合計することができます。このあたりはNumPyの素晴らしいところです。

NumPyを使い偏差平方和を計算

import numpy as np
arr_x = np.array([1, 4, 3, 5, 2])
mean_x = sum(arr_x)/len(arr_x)
sum((arr_x-mean_x)**2) #ブロードキャスト

簡便法による分散の計算

次に簡便法にり計算します。

for文による偏差平方和、分散、標準偏差の計算

はじめに平均を計算し、各データの値から一つ一平均を引き算して偏差を求める必要がないのですっきりします。なお、計算の過程で各データを二乗した値を足しこんでいく必要があるので、square_xという変数を使います。最後にarr_xの合計を2乗してデータ数で割った数値を差し引きます。

for文を使い分散を計算（簡便法）

square_x = 0
for i in list_x:
square_x += i**2
dev_x = square_x-sum(list_x)**2/len(list_x)
var_x = dev_x/len(list_x)
std_x = var_x**0.5
print("xの偏差平方和:{:.1f} 分散:{:.1f} 標準偏差:{:.1f}".format(dev_x, var_x, std_x))

リスト内包表記による方法

リスト内包表記は慣れないと難しいところがありますが、定義通りの方法より簡便法の方がすっきりしていて書きやすくなります。

リスト内包表記を使い偏差平方和を計算（簡便法）

sum([i**2 for i in list_x])-sum(list_x)**2/len(list_x)

ラムダ式とmap関数による方法

ラムダ式も簡便法の方がすっきりしていて書きやすくなります。

ラムダ式～map関数を使い偏差平方和を計算（簡便法）

sum(list(map(lambda i: i**2, list_x)))-sum(list_x)**2/len(list_x)

NumPyモジュールによる方法

NumPyを使うと、リストの2乗を合計するところが、驚くほど簡単になります。

NumPyを使い偏差平方和を計算（簡便法）

import numpy as np
arr_x = np.array([1, 4, 3, 5, 2])
sum(arr_x**2)-sum(arr_x)**2/len(arr_x)

共分散の計算

共分散は2つのデータの関係を計算することから、新たな変数を使います。

用語	変数	意味
yのリスト	list_y	4, 10, 6, 11, 4の5つのデータをリスト形式で定義
積和	product_xy	xとyそれぞれの偏差（平均との差）の2乗の合計
共分散	cov_xy	積和をデータ数nで割った数値

定義通りの分散の計算

for文による積和、共分散の計算

list_xに加えてlist_yを定義して共分散を計算します。2つのリストから同時に1つずつlist_xはi、list_yはjとして取り出して計算するときには、zip関数を使う必要があります。

for文を使い共分散を計算

list_x = list([1, 4, 3, 5, 2])
list_y = list([4, 10, 6, 11, 4])
mean_x = sum(list_x)/len(list_x)
mean_y = sum(list_y)/len(list_y)
product_xy = 0
for i, j in zip(list_x, list_y):
product_xy += (i-mean_x)*(j-mean_y)
cov_xy = product_xy/len(list_x)
print("xyの積和:{:.1f} 共分散:{:.1f}".format(product_xy, cov_xy)) dev_x += (i-mean_x)**2

結果は次の通りです。

xyの積和:20.0  共分散:4.0

リスト内包表記による方法

For文を使う方法と考え方は全く同じです。2つのリストはzip関数で囲う必要があります。

リスト内包表記を使い積和を計算

mean_x = sum(list_x)/len(list_x)
mean_y = sum(list_y)/len(list_y)
sum([(i-mean_x)*(j-mean_y) for i, j in zip(list_x, list_y)])

次のように書いてしまうとlist_xとlist_yの全ての組み合わせ（総当たり）で計算されるので注意が必要です。

[(i-mean_x)*(j-mean_y) for i in list_x for j in list_y]

ラムダ式とmap関数による方法

ラムダ式とmap関数の組み合わせの場合、2つのリストをzip関数で囲う必要はありません。逆に囲ってしまうとエラーになります。

ラムダ式とmap関数を使い積和を計算

mean_x = sum(list_x)/len(list_x)
mean_y = sum(list_y)/len(list_y)
sum((list(map(lambda i, j: (i-mean_x)*(j - mean_y), list_x, list_y)))) # x,y zipなし

NumPyモジュールによる方法

最後にNumPyモジュールを使う方法です。ブロードキャストを使うと本当に簡単です。

NumPyにより積和を計算

import numpy as np
arr_x = np.array([1, 4, 3, 5, 2])
arr_y = np.array([4, 10, 6, 11, 4])
mean_x = np.mean(arr_x)
mean_y = np.mean(arr_y)
product_xy = sum((arr_x-mean_x)*(arr_y-mean_y))
cov_xy = product_xy/len(arr_x)
print("積和:{:.1f} 共分散:{:.1f}".format(product_xy, cov_xy))

簡易法による共分散の計算

簡易法による、分散の計算です。前節と同じように計算してみます。

for文による偏差平方和、分散、標準偏差の計算

計算式通りにプログラミングすれば、問題はないと思います。

for文を使い共分散を計算（簡便法）

square_xy = 0
for i, j in zip(list_x, list_y):
square_xy += i*j
product_xy = square_xy-sum(list_x)*sum(list_y)/len(list_x)
cov_xy = product_xy/len(list_x)
print("xyの積和:{:.1f} 共分散:{:.1f}".format(product_xy, cov_xy))

リスト内包表記による方法

リスト内包表記は慣れないと難しいところがありますが、簡便法の方がすっきりします。

リスト内包表記を使い積和を計算（簡便法）

sum([i*j for i, j in zip(list_x, list_y)])-sum(list_x)*sum(list_y)/len(list_x)

ラムダ式とmap関数による方法

ラムダ式も簡便法の方がすっきりしていて書きやすくなります。

ラムダ式とmap関数を使い積和を計算（簡便法）

sum((list(map(lambda i, j: i*j, list_x, list_y))))-sum(list_x)*sum(list_y)/len(list_x)

NumPyモジュールによる方法

NumPyを使うとリストの2乗を合計するところが、驚くほど簡単になります。

NumPyにより積和を計算（簡便法）

import numpy as np
arr_x = np.array([1, 4, 3, 5, 2])
arr_y = np.array([4, 10, 6, 11, 4])
product_xy = sum((arr_x)*(arr_y))-sum(arr_x)*sum(arr_y)/len(arr_x)
cov_xy = product_xy/len(arr_x)
print("積和:{:.1f} 共分散:{:.1f}".format(product_xy, cov_xy))

NumPyの関数を使い計算する

NumPyには、分散や共分散などを簡単に計算できる関数が用意されています。関数にデータを指定すれば結果を返してくれますが、自由度の部分で注意が必要です。

var関数による分散、std関数による標準偏差の計算

はじめにvar関数は分散、std関数は標準偏差を計算します。これらの関数では配列（ndarrayでもlistでも可）とともにddofを指定する必要があります。計算の対象が母集団であるときは偏差平方和をデータ数で割る方法により”ddof=0”とします。一方、計算の対象が母集団から抜き取ったサンプルであるときは、サンプル数$-$1で割る方法によりddof=1”とします。ここでは、簡単にするためddof=0としました。

NumPyの関数を使って分散を計算する

var_x = np.var(arr_x, ddof=0) # ddof=Delta Degrees of Freedom
std_x = np.std(arr_x, ddof=0)
print("xの分散:{:.1f} 標準偏差:{:.1f}".format(var_x, std_x))
var_y = np.var(arr_y, ddof=0) # ddof=Delta Degrees of Freedom
std_y = np.std(arr_y, ddof=0)
print("yの分散:{:.1f} 標準偏差:{:.1f}".format(var_y, std_y))

結果は次の通りです。

xの分散:2.0  標準偏差:1.4
yの分散:8.8  標準偏差:3.0

cov関数による共分散の計算

共分散はcov関数を使います。考え方はvar関数と同じですが、共分散なので配列を2つ指定します。

np.cov(arr_x, arr_y, ddof=0)

結果は次の通りとなり、初めて見ると少し戸惑います。

array([[2. , 4. ],　　　　
       [4. , 8.8]])

cov関数の結果は次のとおりで表示されます。これは、cov関数が重回帰分析を想定しているためです。

array([[xの分散  , 共分散 ],
       [共分散 ,  yの分散]])

このため、共分散だけを取り出したいときは、次の通りにする必要があります。

np.cov(arr_x, arr_y, ddof=0)[0][1]

arrayの1番目の要素[2,4]の2番目の要素4が共分散になるためです。

Ddofの指定

var関数、std関数、cov関数ではddofを省略してもエラーになりませんが、0か１か明示することをお勧めします。なぜなら、var関数,std関数は何も指定しないときは0ですが、cov関数は1となり、思わぬところでずれが生じてしまうためです。

ちなみに、ddof は”Delta Degrees of Freedom”の略で、統計学でいう自由度をいいます。