Python統計プログラミング
平均・分散・相関係数・グラフ描画（ヒストグラム、散布図）
簡単な統計量とグラフ描画のプログラムをPythonで作ってみた。１変量および２変量の記述統計量、グラフ描画、ファイル入出力の３セッションに分けて説明する。他の統計分析プログラム例は別のサイトに用意している。
なお、グラフにおける日本語表記については、モジュールjapanize_matplotlibによる場合とjapanize_matplotlibを用いない場合について例を上げた。
統計についての解説は、岡本安晴「データを使いこなすための統計入門」を参照されたい。Pythonのインストールについては、このウェブサイトで説明している。
Pythonによるデータ分析の解説書として
岡本安晴（著）「いまさら聞けないPythonでデータ分析―多変量解析，ベイズ統計分析（PyStan，PyMC）―」丸善出版
を用意した。
なお、本ウェブサイトのプログラム例は圧縮ファイルSmplStatPy.zipにまとめたが、ダウンロード・展開して、自由に御利用下さい。
記述統計量
１変量の場合
表１の５０個の値の平均値と分散などを求めてみる。
|
表１ ５０個のスコア |
|
66, 59, 62, 64, 63,68, 65, 59, 68, 64, 65, 51, 67, 64, 83,59, 61, 62, 57, 72, 65, 64, 54, 60, 53,65, 67, 60, 53, 79, 74, 53, 61, 68, 75,50, 57, 55, 66, 56, 55, 61, 70, 71, 49,69, 70, 80, 73, 72 |
プログラムStatOneVar.pyを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝ StatOneVar.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
def sqr(x):
return x * x
data = [66, 59, 62, 64,
63,
68,
65, 59, 68, 64,
65,
51, 67, 64, 83,
59,
61, 62, 57, 72,
65,
64, 54, 60, 53,
65,
67, 60, 53, 79,
74,
53, 61, 68, 75,
50,
57, 55, 66, 56,
55,
61, 70, 71, 49,
69,
70, 80, 73, 72]
n = len(data) # The number of scores
print("n
=%d"%n)
print("Data...")
i = 1
for x in data:
print("%5d:"%i +
" %f"%x)
i+=1
sum = 0.0
for x in data:
sum += x
mean = sum / n
print("Mean =
%f"%mean)
ssum = 0.0
for x in data:
ssum += sqr(x - mean)
var = ssum / n
sd = math.sqrt(var)
print("Variance
= %f"%var + "
sd = %f"%sd)
uvar = ssum / (n -
1.0)
sd_uvar =
math.sqrt(uvar)
print("Unbiased
Var. = %f"%uvar + " sd(uvar) =
%f"%sd_uvar)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
表１のデータをリストdataとして設定（assignment記号 = の働きについては、本節末に補足説明を用意した。）している。平均値、分散、不偏分散、および分散あるいは不偏分散の平方根としての標準偏差が算出されている。それぞれは、次式
、
、
、
、
で算出される。
プログラムStatOneVar.pyを実行すると、図１のような結果になる。
図１
２変量の場合
体重と身長のように２つの変量がペアになっている場合について考える。
|
表２ ［体重, 身長］の５０人分のデータ |
|
[70, 168], [68, 165], [62, 159], [56, 156], [70, 170], [64, 162], [59, 162], [70, 169], [68, 168], [69, 170], [67, 168], [63, 166], [58, 163], [74, 170], [65, 168], [60, 160], [65, 165], [72, 175], [61, 164], [77, 174], [68, 164], [64, 161], [70, 166], [47, 147], [66, 165], [68, 170], [73, 174], [74, 172], [63, 164], [83, 182], [58, 159], [64, 166], [67, 167], [75, 174], [67, 163], [68, 166], [68, 167], [66, 165], [63, 165], [55, 154], [61, 163], [60, 160], [62, 164], [62, 160], [59, 160], [51, 155], [59, 163], [68, 170], [56, 157], [79, 180] |
プログラムStatTwoVars.pyでは、表２のデータがリストpairsとして設定されている。
＝＝＝＝＝＝＝＝＝＝＝＝＝ StatTwoVars.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
def sqr(x):
return x * x
pairs = [[70, 168], [68,
165], [62, 159], [56, 156], [70, 170],
[64, 162], [59, 162], [70, 169], [68, 168], [69, 170],
[67, 168], [63, 166], [58, 163], [74, 170], [65, 168],
[60, 160], [65, 165], [72, 175], [61, 164], [77, 174],
[68, 164], [64, 161], [70, 166], [47, 147], [66, 165],
[68, 170], [73, 174], [74, 172], [63, 164], [83, 182],
[58, 159], [64, 166], [67, 167], [75, 174], [67, 163],
[68, 166], [68, 167], [66, 165], [63, 165], [55, 154],
[61, 163], [60, 160], [62, 164], [62, 160], [59, 160],
[51, 155], [59, 163], [68, 170], [56, 157], [79, 180]]
n = len(pairs)
print("n =
%d"%n)
sum_x = 0.0
sum_y = 0.0
for i in range(0, n):
sum_x += pairs[i][0]
sum_y += pairs[i][1]
mean_x = sum_x / n
mean_y = sum_y / n
print("Mean_x =
%f"%mean_x)
print("Mean_y =
%f"%mean_y)
ssum_x = 0.0
ssum_y = 0.0
sum_xy = 0.0
for i in range(0, n):
ssum_x += sqr(pairs[i][0] -
mean_x)
ssum_y += sqr(pairs[i][1] -
mean_y)
sum_xy += (pairs[i][0] - mean_x)
* (pairs[i][1] - mean_y)
var_x = ssum_x / n
var_y = ssum_y / n
print("Var_x =
%f"%var_x)
print("Var_y =
%f"%var_y)
sd_x =
math.sqrt(var_x)
sd_y =
math.sqrt(var_y)
print("sd_x =
%f"%sd_x)
print("sd_y =
%f"%sd_y)
cor = sum_xy /
math.sqrt(ssum_x * ssum_y)
print("r =
%f"%cor)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
２変量の関係を表すものとして相関係数（correlation coefficient）が求められている。これは、次式
によって算出されている。
プログラムStatTwoVars.pyを実行すると図２のような結果になる。
図２
図２の結果は、Windowsのコンソール（端末）での実行結果であるが、図１のようにPythonのIDLE上で実行することもできる。
Assignment = については、以下のコードに示されるような注意が必要である。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
a = [1, 2]
b = a
# Check-11
b[0] = 3
# Check-12
print("Check-1")
print(a)
print(b)
a1 = [1, 2]
b1 = a1
# Check-21
a1[0] = 5
# Check-22
print("Check-2")
print(b1)
print(a1)
a2 = [1, 2]
b2 = a2
# Check-31
b2 = [6, 7]
# Check-32
print("Check-3")
print(a2)
print(b2)
a3 = 10
b3 = a3
b3 = 15
print("Check-4")
print(a3)
print(b3)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のコードを実行すると以下の結果を得る。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Check-1
[3, 2]
[3, 2]
Check-2
[5, 2]
[5, 2]
Check-3
[1, 2]
[6, 7]
Check-4
10
15
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
「Check-1」に続く出力では、同じ出力「[3, 2]」が得られている。これは、Check-11における代入文によって、変数bが変数aと同じものを指し示すものとなっていることを意味している。
「Check-2」に続く出力も、代入がb1ではなくa1（b1 = a1の右辺）の要素に行われているが、2つの変数b1とa1が同じものを表していることを意味していることは変数aとbの場合と同様であることを意味するものである。
「Check-3」に続く出力では、変数b2の値の変更は変数a2には及んでいない。これは、変数の値全体を変更する場合は、変数は新しい変更値を指し示すものとして設定され、それまでの同じものを指し示していた変数a2から切り離されたことを意味している。
「Check-4」に続く出力は整数の場合を扱っているが、整数はそれ１つで変数の表す値全体であるので、変数b3に値を設定すると、それまで同じものを表していた変数a3とは別のものを指し示すことになる。
グラフ描画
上の表１のデータに対して、ヒストグラムを描くプログラムHistogram_1.pyを作成してみた。描画プログラムではパッケージmatplotlibを用いている（importしている）が、これはPythonに追加インストールするものである。追加インストールの方法は＜別のウェブサイト＞で説明している。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Histogram_1.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import
matplotlib.pyplot as mpt
import math
def sqr(x):
return x * x
data = [66, 59, 62,
64, 63,
68,
65, 59, 68, 64,
65,
51, 67, 64, 83,
59,
61, 62, 57, 72,
65,
64, 54, 60, 53,
65,
67, 60, 53, 79,
74,
53, 61, 68, 75,
50,
57, 55, 66, 56,
55,
61, 70, 71, 49,
69,
70, 80, 73, 72]
n = len(data)
# The number of scores
print("n
=%d"%n)
print("Data...")
i = 1
for x in data:
print("%5d:"%i +
" %f"%x)
i+=1
mpt.title("Histogram") # The title of the histogram
mpt.xlabel("Score")
# The label of the
x-axis
mpt.ylabel("Frequency") # The label of the y-axis
mpt.hist(data, bins =
10) # Making the histogram
mpt.show()
# Show the histogram
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
モジュールmatplotlibを用いているが、これはコンソール（図２に示されているようなWindow）において
pip install matplotlib
とコマンドpipを実行してインストールした。
ヒストグラムの描画は、最後の次のコード
mpt.title("Histogram") # The title of the histogram
mpt.xlabel("Score")
# The label of the
x-axis
mpt.ylabel("Frequency") # The label of the y-axis
mpt.hist(data,
bins = 10) # Making the histogram
mpt.show()
# Show the histogram
によって行っている。show()の実行によって図３のヒストグラムが表示される。なお、グラフの日本語表記については、本ウェブサイトの後半で説明している。
図３
下のフロッピーディスクのアイコンのボタンをクリックすると、表示画像を保存することができる。右上のＸ印のクリックで、表示画面は閉じられ、show()の実行終了となる。
図３のヒストグラムは、枠との関係を調整した方が良さそうである。この調整機能を付け加えたプログラムHistogram_2.pyを以下のように作成してみた、
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Histogram_2.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import
matplotlib.pyplot as mpt
import math
def sqr(x):
return x * x
data = [66, 59, 62,
64, 63,
68,
65, 59, 68, 64,
65,
51, 67, 64, 83,
59,
61, 62, 57, 72,
65,
64, 54, 60, 53,
65,
67, 60, 53, 79,
74,
53, 61, 68, 75,
50,
57, 55, 66, 56,
55,
61, 70, 71, 49,
69,
70, 80, 73, 72]
n = len(data)
# The number of scores
print("n
=%d"%n)
print("Data...")
i = 1
for x in data:
print("%5d:"%i +
" %f"%x)
i+=1
sum = 0.0
for x in data:
sum += x
mean = sum /
len(data)
print("Mean =
%f"%mean)
ssum = 0.0
for x in data:
ssum += sqr(x - mean)
var = ssum /
len(data)
sd = math.sqrt(var)
print("Variance
= %f"%var + "
sd = %f"%sd)
uvar = ssum /
(len(data) - 1.0)
sd_uvar =
math.sqrt(uvar)
print("Unbiased
Var. = %f"%uvar + " sd(uvar) =
%f"%sd_uvar)
#
# Check the
minimum and maximum values
#
vmin = data[0]
n_max = n
vmax=data[0]
for x in data:
if vmin > x: vmin = x
if vmax < x: vmax = x
print("min =
%f"%vmin + " max =
%f"%vmax)
n_bin = 10
while(True):
mpt.title("Histogram")
# The title of the histogram
mpt.xlabel("Score")
# The label of the
x-axis
mpt.ylabel("Frequency")
# The label of the y-axis
mpt.axis([vmin,vmax,0,n_max])
# Set the boundaries
of the axes
mpt.hist(data, bins =
n_bin) # Make the histogram
mpt.show()
# Show the histogram
n_bin =
int(input("n_bin/%d = "%n_bin)) # Set the number of
if (n_bin < 2):
break
# The program ended
n_max =
int(input("n_max/%d = "%n_max)) # The upper boundary of the y-axis
vmin =
float(input("x_min/%f = "%vmin)) # The left boundary of the x-axis
vmax =
float(input("x_max/%f = "%vmax)) # The right boundary of the x-axis
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
以下のwhile文内の繰り返し処理において枠の調整を行っている。
while(True):
mpt.title("Histogram")
# The title of the
histogram
mpt.xlabel("Score")
# The label of the
x-axis
mpt.ylabel("Frequency") # The label of the y-axis
mpt.axis([vmin,vmax,0,n_max])
# Set the boundaries
of the axes
mpt.hist(data, bins =
n_bin) # Make the histogram
mpt.show()
# Show the histogram
n_bin =
int(input("n_bin/%d = "%n_bin)) # Set the number of
if (n_bin < 2):
break
# The program ended
n_max =
int(input("n_max/%d = "%n_max)) # The upper boundary of the y-axis
vmin =
float(input("x_min/%f = "%vmin)) # The left boundary of the x-axis
vmax =
float(input("x_max/%f = "%vmax)) # The right boundary of the x-axis
上のプログラムHistogram_2.pyを実行すると、図４のヒストグラムが表示される。
図４
画面の右上のＸ印をクリックして画面を閉じると、Python IDLEのコンソールが入力待ちの状態になる（図５）。
図５
適宜、適当な数値を入力していくが、ビンの数（ヒストグラムの区間数）として１以下の数を設定すると、プログラムの実行終了となる。
図６
直前の描画において用いられた値は、／の右側に表示されている。これらを参考にして、適当な値を入力していく（図６）。図６の状態でEnterキーを押すと、図７のヒストグラムが表示される。
図７
２変量の場合
散布図を描くプログラムScatter.pyを以下のように作成してみた。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Scatter.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
import
matplotlib.pyplot as plt
def sqr(x):
return x * x
#
# Data
#
pairs = [[70, 168], [68,
165], [62, 159], [56, 156], [70, 170],
[64, 162], [59, 162], [70, 169], [68, 168], [69, 170],
[67, 168], [63, 166], [58, 163], [74, 170], [65, 168],
[60, 160], [65, 165], [72, 175], [61, 164], [77, 174],
[68, 164], [64, 161], [70, 166], [47, 147], [66, 165],
[68, 170], [73, 174], [74, 172], [63, 164], [83, 182],
[58, 159], [64, 166], [67, 167], [75, 174], [67, 163],
[68, 166], [68, 167], [66, 165], [63, 165], [55, 154],
[61, 163], [60, 160], [62, 164], [62, 160], [59, 160],
[51, 155], [59, 163], [68, 170], [56, 157], [79, 180]]
n = len(pairs) # The number
of pairs (xi, yi)
print("n =
%d"%n)
sum_x = 0.0
sum_y = 0.0
for i in range(0, n):
sum_x += pairs[i][0]
sum_y += pairs[i][1]
mean_x = sum_x / n
mean_y = sum_y / n
print("Mean_x =
%f"%mean_x)
print("Mean_y =
%f"%mean_y)
ssum_x = 0.0
ssum_y = 0.0
sum_xy = 0.0
for i in range(0, n):
ssum_x += sqr(pairs[i][0] -
mean_x)
ssum_y += sqr(pairs[i][1] - mean_y)
sum_xy += (pairs[i][0] -
mean_x) * (pairs[i][1] - mean_y)
var_x = ssum_x /
n
# Variance of X
var_y = ssum_y /
n
# Variance of Y
print("Var_x =
%f"%var_x)
print("Var_y =
%f"%var_y)
sd_x = math.sqrt(var_x)
# Standard deviation
of X
sd_y =
math.sqrt(var_y)
# Standard deviation
of Y
print("sd_x =
%f"%sd_x)
print("sd_y =
%f"%sd_y)
cor = sum_xy /
math.sqrt(ssum_x * ssum_y) # The correlation coefficient
print("r =
%f"%cor)
data_x = [] # Values of the coordinate Xs of the
pairs
data_y = [] # Values of the coordinate Ys of the
pairs
for i in range(0, n):
data_x.append(pairs[i][0])
data_y.append(pairs[i][1])
print(data_x)
print(data_y)
plt.title("Scattergram...r
= %6.3f"%cor)
# The title of the
scattergram
plt.xlabel("Weight")
# The label of the
x-axis
plt.ylabel("Height")
# The label of the
y-axis
plt.plot(data_x,
data_y, 'bo')
# Make the
scattergram. 'bo' -> blue circle
plt.show()
# Show the scattergram
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
散布図描画の関数plotは、以下のように、X座標のリストdata_xとＹ座標のリストdata_yに分けて別個のパラメータとして渡している。
plt.plot(data_x, data_y, 'bo') # Make the scattergram. 'bo' -> blue circle
第３パラメータ「'bo’」は、散布図の点を青の円で描くことを指示するものである。
Ｘ座標のリストdata_xとＹ座標のリストdata_yを、［Ｘ座標、Ｙ座標］のペアでデータが設定されているリストpairsから作成するために以下のコードを用いている。
data_x = [] # Values of the coordinate Xs of the pairs
data_y = [] # Values of the coordinate Ys of the pairs
for i in range(0, n):
data_x.append(pairs[i][0])
data_y.append(pairs[i][1])
プログラムScatter.pyを実行すると図８の散布図が表示される。なお、グラフの日本語表記については、本ウェブサイトの後半で説明している。
図８
ファイル入出力
入力データは、プログラム中に書き込むよりは別ファイルとして用意する方が便利である。また、出力もコンソールではなく、ファイルに出力した方が後で使いやすい。
１変量データの場合
表１のデータを、DataOneVar.txtとして以下のようなテキストファイルDataOneVar.txtとして用意した。このファイルは圧縮ファイルSmplStatPy.zipに含まれている。
＝＝＝＝＝＝＝＝＝＝＝＝ DataOneVar.txt ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
/
66
59
62
64
63
68
・
・
・
49
69
70
80
73
72
/
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
ファイル入出力に対応したプログラムStatOneVarF.pyは以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ StatOneVarF.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
import
matplotlib.pyplot as mpt
def sqr(x):
return x * x
#
# Prepare
the input data file
#
s_in =
input("Input data file = ")
f_in = open(s_in,
"r")
#
# Set the contents
of the input data file in the object data
#
data_f =
f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
s_out =
input("Output data file = ")
f_out = open(s_out,
"w")
f_out.write("Input
data file = " + s_in + "\n\n")
pos = 0
while True:
if
data_f[pos].find("/", 0) == 0: break
pos += 1
n = 0
data = []
while True:
pos += 1
if
data_f[pos].find("/", 0) == 0: break
data.append(float(data_f[pos]))
n += 1
print("Data...")
i = 1
for x in data:
f_out.write("%5d:"%i
+ " %f\n"%x)
i+=1
f_out.write("\n")
sum = 0.0
for x in data:
sum += x
mean = sum / n
print("Mean =
%f"%mean)
f_out.write("Mean
= %f\n"%mean)
ssum = 0.0
for x in data:
ssum += sqr(x - mean)
var = ssum / n
sd = math.sqrt(var)
print("Variance
= %f"%var + "
sd = %f"%sd)
f_out.write("Variance
= %f"%var + "
sd = %f\n"%sd)
uvar = ssum / (n -
1.0)
sd_uvar =
math.sqrt(uvar)
print("Unbiased
Var. = %f"%uvar + " sd(uvar) =
%f"%sd_uvar)
f_out.write("Unbiased
Var. = %f"%uvar + " sd(uvar) =
%f\n"%sd_uvar)
f_out.close()
print(s_out + "
was saved.")
vmin = data[0]
n_max = n
vmax=data[0]
for x in data:
if vmin > x: vmin = x
if vmax < x: vmax = x
print("min =
%f"%vmin + " max =
%f"%vmax)
n_bin = 10
while(True):
mpt.title("Histogram")
mpt.xlabel("Scores")
mpt.ylabel("Frequencies")
mpt.axis([vmin,vmax,0,n_max])
mpt.hist(data, bins = n_bin)
mpt.show()
n_bin =
int(input("n_bin/%d = "%n_bin))
if (n_bin < 2): break
n_max =
int(input("n_max/%d = "%n_max))
vmin =
float(input("x_min/%f = "%vmin))
vmax =
float(input("x_max/%f = "%vmax))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
次のコードによって、コンソールから入力した名前のファイルを入力ファイルとしている。
s_in = input("Input data file = ")
f_in = open(s_in, "r")
出力ファイルは以下のコードによって設定している。
s_out = input("Output data file = ")
f_out = open(s_out, "w")
入力データファイルの内容は、コード
data_f = f_in.readlines()
によってリストdata_fに設定されている。入力データファイルの構成（これはリストdata_fに文字列のリストとして取り込まれている）は、データの書かれている部分がスラッシュ / で始まる行で挟まれているので、以下のコードによってデータ（data_fの内容）をリストdataに設定している。
pos = 0
while True:
if data_f[pos].find("/", 0) == 0: break
pos += 1
n = 0
data = []
while True:
pos += 1
if data_f[pos].find("/", 0) == 0: break
data.append(float(data_f[pos]))
n += 1
この後の処理は、上のプログラムHistogram_2.pyに対応している。
プログラムを実行すると、まず、入力データファイル名を聞いてくる（図９）。
図９
入力データファイル名の設定後、Enterキーを押すと、次に出力ファイル名を聞いてくる（図１０）。
図１０
出力ファイル名として適当なテキストファイル名を設定してEnterキーを押すと計算が始まり、図１１の画面が表示される。
図１１
これは図４と同じである。操作も図４に対するものと同じように行えばよい。
プログラムの実行終了後、出力ファイルを開くと、入力データファイルDataOneStat.txtの場合は、出力ファイルは以下のようになっている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ResultsOneVar.txt ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Input data file =
DataOneVar.txt
1: 66.000000
2: 59.000000
3: 62.000000
4: 64.000000
5: 63.000000
・
・
・
47: 70.000000
48: 80.000000
49: 73.000000
50: 72.000000
Mean = 63.680000
Variance =
60.217600 sd =
7.760000
Unbiased Var. =
61.446531
sd(uvar) = 7.838784
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
２変量データの場合
表２のデータから、以下のようにテキストファイルDataTwoVars.txtを作成した。このファイルも圧縮ファイルSmplStatPy.zipに含まれている。
＝＝＝＝＝＝＝＝＝＝＝＝＝ DataTwoVars.txt ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
/
70 168
68 165
62 159
56 156
70 170
64 162
59 162
・
・
・
62 160
59 160
51 155
59 163
68 170
56 157
79 180
/
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
データがスラッシュ / で始まる行で挟まれているので、読み込み方は１変量のプログラムStatOneVarF.pyの場合と基本的に同じである。２変量用のプログラムScatterF.pyを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ScatterF.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
import
matplotlib.pyplot as mpt
def sqr(x):
return x * x
#
# Prepare
the input data file
#
s_in =
input("Input data file = ")
f_in = open(s_in,
"r")
#
# Set the
contents of the input data file in the object data
#
data_f =
f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
s_out =
input("Output data file = ")
f_out = open(s_out,
"w")
f_out.write("Input
data file = " + s_in + "\n\n")
pos = 0
while True:
if data_f[pos].find("/", 0) ==
0: break # Jump to teh first line with slash
at the head
pos += 1
n = 0
pairs = []
while True:
pos += 1
if
data_f[pos].find("/", 0) == 0: break # Exit when the line with slash at the
head is met
#
# Add the
data
#
temp_s = data_f[pos].split()
pairs.append([])
pairs[n].append(float(temp_s[0]))
pairs[n].append(float(temp_s[1]))
n += 1
f_out.write("\n")
f_out.write("The
number of pairs = %d\n"%n)
f_out.write("Data...\n")
i = 1
for v in pairs:
f_out.write("%5d: %f
%f\n"%(i, v[0], v[1]))
i += 1
sum_x = 0.0
sum_y = 0.0
for i in range(0, n):
sum_x += pairs[i][0]
sum_y += pairs[i][1]
mean_x = sum_x / n
mean_y = sum_y / n
print("Mean_x =
%f"%mean_x)
print("Mean_y =
%f"%mean_y)
f_out.write("\nMean_x
= %f\n"%mean_x)
f_out.write("Mean_y
= %f\n"%mean_y)
ssum_x = 0.0
ssum_y = 0.0
sum_xy = 0.0
for i in range(0, n):
ssum_x += sqr(pairs[i][0] -
mean_x)
ssum_y += sqr(pairs[i][1] -
mean_y)
sum_xy += (pairs[i][0] - mean_x) *
(pairs[i][1] - mean_y)
var_x = ssum_x /
n
# Variance of X
var_y = ssum_y /
n
# Variance of Y
print("Var_x =
%f"%var_x)
print("Var_y =
%f"%var_y)
f_out.write("Var_x
= %f\n"%var_x)
f_out.write("Var_y
= %f\n"%var_y)
sd_x =
math.sqrt(var_x)
# Standard deviation
of X
sd_y =
math.sqrt(var_y)
# Standard deviation
of Y
print("sd_x =
%f"%sd_x)
print("sd_y =
%f"%sd_y)
f_out.write("sd_x
= %f\n"%sd_x)
f_out.write("sd_y
= %f\n"%sd_y)
cor = sum_xy /
math.sqrt(ssum_x * ssum_y)
# The correlation
coefficient
print("r =
%f"%cor)
f_out.write("r =
%f\n"%cor)
data_x = [] # Values of the coordinate Xs of the
pairs
data_y = [] # Values of the coordinate Ys of the
pairs
for i in range(0, n):
data_x.append(pairs[i][0])
data_y.append(pairs[i][1])
print(data_x)
print(data_y)
f_out.close()
print("Output
file " + s_out + " was saved.")
mpt.title("Scattergram...r
= %6.3f"%cor)
mpt.xlabel("Weight")
mpt.ylabel("Height")
mpt.plot(data_x,
data_y, 'g^')
# 'g^' -> green
triangle
mpt.show()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
１行に１組のデータ、ＸとＹ、が含まれているので、それを分離してリストpairsに格納するため、以下のコードが用いられている。
temp_s = data_f[pos].split()
pairs.append([])
pairs[n].append(float(temp_s[0]))
pairs[n].append(float(temp_s[1]))
また、上のプログラムScatterF.pyでは、散布図の点が緑色の三角形で描かれるように、第３パラメータの値が「'g^’」となっている。
mpt.plot(data_x, data_y, 'g^') # 'g^' -> green triangle
このとき、散布図は図１２のように描かれる。
図１２
出力ファイを図１３の名前ResultsTwoVars.txtで設定すると、
図１３
その内容は以下に示すようになっている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ResultsTwoVars.txt ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Input data file =
DataTwoVars.txt
The number of pairs =
50
Data...
1: 70.000000 168.000000
2: 68.000000 165.000000
3: 62.000000 159.000000
4: 56.000000 156.000000
5: 70.000000 170.000000
・
・
・
45: 59.000000 160.000000
46: 51.000000 155.000000
47: 59.000000 163.000000
48: 68.000000 170.000000
49: 56.000000 157.000000
50: 79.000000 180.000000
Mean_x = 65.240000
Mean_y = 165.300000
Var_x = 46.822400
Var_y = 40.130000
sd_x = 6.842690
sd_y = 6.334824
r = 0.942805
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上の計算結果は、図２に示されている計算結果と同じであることが確認できる。
グラフの日本語表示（モジュールjapanize_matplotlibによる場合）
グラフ描画のモジュールmatplotlibに加えて、モジュールjapanize_matplotlibをimportすると日本語が表示できる。
モジュールjapanize_matplotlibはコマンドpipでインストールできる。
pip install
japanize_matplotlib
と実行すればよい。しかし、Anacondaのcondaコマンドではインストールできなかった（2020.0613）が、pipコマンドでインストールすることができた。
ヒストグラムの日本語表記のサンプルプログラムは以下のようである。
import matplotlib.pyplot as mpt
import japanize_matplotlib
import math
def sqr(x):
return x * x
data = [66, 59, 62, 64, 63,
68,
65, 59, 68, 64,
65,
51, 67, 64, 83,
59,
61, 62, 57, 72,
65,
64, 54, 60, 53,
65,
67, 60, 53, 79,
74,
53, 61, 68, 75,
50,
57, 55, 66, 56,
55,
61, 70, 71, 49,
69,
70, 80, 73, 72]
n = len(data)
# The number of scores
print("個数 = {}".format(n))
print("データ...")
i = 1
for x in data:
print("{0:}: {1:}".format(i, x))
i+=1
mpt.title("ヒストグラム") # タイトル
mpt.xlabel("スコア")
# X軸のラベル
mpt.ylabel("度数")
# Y軸のラベル
mpt.hist(data, bins = 10) # ヒストグラムの描画
mpt.show()
# 描画グラフの表示
上のスクリプトを実行すると、下図のグラフが表示される。
散布図の場合のプログラム例を以下に示す。
import math
import matplotlib.pyplot as plt
import japanize_matplotlib
def sqr(x):
return x * x
#
# Data
#
pairs = [[70, 168], [68, 165], [62,
159], [56, 156], [70, 170],
[64, 162], [59, 162], [70, 169], [68, 168], [69, 170],
[67, 168], [63, 166], [58, 163], [74, 170], [65, 168],
[60, 160], [65, 165], [72, 175], [61, 164], [77, 174],
[68, 164], [64, 161], [70, 166], [47, 147], [66, 165],
[68, 170], [73, 174], [74, 172], [63, 164], [83, 182],
[58, 159], [64, 166], [67, 167], [75, 174], [67, 163],
[68, 166], [68, 167], [66, 165], [63, 165], [55, 154],
[61, 163], [60, 160], [62, 164], [62, 160], [59, 160],
[51, 155], [59, 163], [68, 170], [56, 157], [79, 180]]
n = len(pairs) # データ(Xi, Yi)の個数
print("個数 = {}".format(n))
sum_x = 0.0
sum_y = 0.0
for i in range(0, n):
sum_x += pairs[i][0]
sum_y += pairs[i][1]
mean_x = sum_x / n
mean_y = sum_y / n
print("Xの平均 =
{0:.2f}".format(mean_x))
print("Yの平均 =
{0:.2f}".format(mean_y))
ssum_x = 0.0
ssum_y = 0.0
sum_xy = 0.0
for i in range(0, n):
ssum_x += sqr(pairs[i][0] -
mean_x)
ssum_y += sqr(pairs[i][1] -
mean_y)
sum_xy += (pairs[i][0] -
mean_x) * (pairs[i][1] - mean_y)
var_x = ssum_x / n
# Xの分散
var_y = ssum_y / n
# Yの分散
print("Xの分散 =
{0:.2f}".format(var_x))
print("Yの分散 =
{0:.2f}".format(var_y))
sd_x = math.sqrt(var_x)
# Xの標準偏差
sd_y = math.sqrt(var_y)
# Yの標準偏差
print("Xの標準偏差 =
{0:.2f}".format(sd_x))
print("Yの標準偏差 =
{0:.2f}".format(sd_y))
cor = sum_xy / math.sqrt(ssum_x *
ssum_y)
# 相関係数
print("相関係数 = {0:.3}".format(cor))
data_x = [] # データ(X, Y)のXを設定
data_y = [] # データ(X< Y)のYを設定
for i in range(0, n):
data_x.append(pairs[i][0])
data_y.append(pairs[i][1])
print(data_x)
print(data_y)
plt.title("散布図...r =
{0:.3}".format(cor), fontsize = 16) # グラフのタイトル
plt.xlabel("体重", fontsize = 14) # X軸のラベル
plt.ylabel("身長", fontsize = 14) # Y軸のラベル
plt.plot(data_x, data_y, 'bo')
# 散布図の作成. 'bo' -> 青の小円
plt.show()
# グラフの表示
上のプログラムを実行すると下図のグラフが表示される。
グラフの日本語表示（モジュールjapanize_matplotlibによらない場合）
モジュールjapanize_matplotlibを用いず、matplotlib.pyplotの関数rcParamsによってフォントの設定を行う例を以下に示す。
import scipy.stats as ss
import matplotlib.pyplot as plt
d_r0 = ss.multivariate_normal.rvs([0,
0], [[1, 0], [0, 1]], size = 100)
d_r09 =
ss.multivariate_normal.rvs([0,0], [[1, 0.9], [0.9, 1]], size = 100)
plt.rcParams['font.family'] = 'Yu
Gothic'
plt.figure(figsize = (7, 7))
plt.title('散布図', fontsize = 30)
plt.plot(d_r0[:,0], d_r0[:, 1], 'go',
label = '相関係数 = 0')
plt.plot(d_r09[:, 0], d_r09[:, 1],
'bo', label = '相関係数
= 0.9')
plt.xlabel('X軸', fontsize = 20)
plt.ylabel('Y軸', fontsize = 20)
plt.xlim(-3.5, 3.5)
plt.ylim(-3.5, 3.5)
plt.legend(loc = 'lower right', fontsize
= 16)
plt.show()
上の例では、Windowsでの場合なので
plt.rcParams['font.family']
= 'Yu Gothic'
としているが、他のOS、例えばubuntuであれば、他のフォントを指定して
plt.rcParams['font.family']
= 'Noto Sans CJK JP'
などとすればよい。
複数のOSで共通に使うときは、可能性のあるフォントを並べて、
plt.rcParams['font.family']
= ['Yu Gothic', ‘Noto
Sans CJK JP’, ‘Hiragino Maru
Gothic Pro’]
などとすればよい。ただし、この方法は、インストールされていないフォントに対しては、その旨のメッセージが出力されるようになった（2023.02.23）。
インストールされているフォントを調べる方法については、このウェブサイトを参照されたい。
上のリストを実行すると下図のグラフが表示される。
