GCIデータサイエンティスト育成講座
「GCIデータサイエンティスト育成講座」は、東京大学(松尾研究室)が開講している"実践型のデータサイエンティスト育成講座およびDeep Learning講座"で、演習パートのコンテンツがJupyterNoteBook形式で公開(CC-BY-NC-ND)されています。
Chapter2は「Numpy、Scipy、Pandas、Matplotlibの基礎」で、科学分野で多く利用されている数値処理および可視化(グラフ化)ライブラリの使い方を学習していきます。
日本語で学べる貴重で素晴らしい教材を公開いただいていることへの「いいね！」ボタンの代わりに、解いてみた解答を載せてみます。間違っているところがあったらご指摘ください。
Chapter2 Numpy、Scipy、Pandas、Matplotlibの基礎
2.1.1
[やってみよう]
seed(0)の0を変えたり、ランダム抽出の数を増やしたりして、結果どうなっているか表示してみましょう。
random.seed(1)
# 正規分布（平均0、分散1）の乱数を10個発生
print("乱数10個の配列1-1:", random.randn(10))
print("乱数10個の配列1-2:", random.randn(10))
random.seed(0)
# 正規分布（平均0、分散1）の乱数を10個発生
print("乱数10個の配列2-1:", random.randn(10))
print("乱数10個の配列2-2:", random.randn(10))
random.seed(1)
# 元の設定に戻す
> 乱数10個の配列1-1: [ 1.624 -0.612 -0.528 -1.073 0.865 -2.302 1.745 -0.761 0.319 -0.249]
> 乱数10個の配列1-2: [ 1.462 -2.06 -0.322 -0.384 1.134 -1.1 -0.172 -0.878 0.042 0.583]
> 乱数10個の配列2-1: [ 1.764 0.4 0.979 2.241 1.868 -0.977 0.95 -0.151 -0.103 0.411]
> 乱数10個の配列2-2: [ 0.144 1.454 0.761 0.122 0.444 0.334 1.494 -0.205 0.313 -0.854]
<練習問題 1>
1から50までの自然数の和を計算するプログラムを書いて、最後の計算結果を表示させるプログラムを書いてください。ただし、Numpyを使ってください。
exam11Array = np.arange(50) + 1
print(np.sum(exam11Array))
> 1275
<練習問題 2>
正規分布に従う乱数を10個発生させて配列を作成してください。また、その中での最小値、最大値、合計を求めるプログラムを書いてください。
import numpy.random as random
random.seed(1)
exam12Array = random.randn(10)
print(exam12Array.min())
print(exam12Array.max())
print(exam12Array.sum())
> -2.30153869688
> 1.74481176422
> -0.971408908061
<練習問題 3>
要素がすべて3の5行5列の行列を作成し、その行列の2乗をする計算をしてみましょう。
exam13Array = np.full([5,5],3)
print(exam13Array)
print(exam13Array**2) #2乗ってどっちだろう？
print(np.dot(exam13Array, exam13Array)) #2乗ってどっちだろう？
> [[3 3 3 3 3]
> [3 3 3 3 3]
> [3 3 3 3 3]
> [3 3 3 3 3]
> [3 3 3 3 3]]
> [[9 9 9 9 9]
> [9 9 9 9 9]
> [9 9 9 9 9]
> [9 9 9 9 9]
> [9 9 9 9 9]]
> [[45 45 45 45 45]
> [45 45 45 45 45]
> [45 45 45 45 45]
> [45 45 45 45 45]
> [45 45 45 45 45]]
2.1.2
[やってみよう]
上の関数を変更したりして、最小値等を計算を実行してみましょう。
# 関数の定義
def sample_function(x):
return (x**2 + 2*x + 1 + np.sin(x))
# 計算実行
print(minimize_scalar(sample_function,method="Brent"))
> fun: -0.89228084521426954
> nfev: 9
> nit: 8
> success: True
> x: -1.1871514382985477
<練習問題 1>
以下の行列について、行列式を求めてください。
A = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 3 & 2 \\ 3 & 1 & 2 \end{array} \right)
exam21Array = np.array([[1,2,3],[1,3,2],[3,1,2]])
print(linalg.det(exam21Array))
> -12.000000000000002
<練習問題 2>
上と同じ行列について、逆行列、固有値と固有ベクトルを求めてください。
print(linalg.inv(exam21Array))
eig_value, eig_vector = linalg.eig(sample_matrix_data)
print(eig_value)
print(eig_vector)
> [[-0.333 0.083 0.417]
> [-0.333 0.583 -0.083]
> [ 0.667 -0.417 -0.083]]
> [-1.+0.j 2.+0.j 2.+0.j]
> [[ 0.577 -0.816 0.428]
> [ 0.577 0.408 -0.816]
> [ 0.577 0.408 0.389]]
<練習問題 3>
以下の関数が0となる解を求めてみましょう。
\begin{eqnarray} f(x) = x^3 + 2x+ 1 \end{eqnarray}
def examFunc(x):
return x ** 3 + 2 * x + 1
print(newton(examFunc,0))
> -0.4533976515164037
2.1.3
[やってみよう]
他にも条件を変更（birth_yearが1990未満など）して、フィルターを実行してみましょう。
attri_data_frame1[attri_data_frame1['birth_year'] < 1990]
[やってみよう]
他にも変数を変えて、実行してみましょう。Englishの場合はどうなりますか。
attri_data_frame2.groupby("sex")["English"].mean()
> sex
> F 25.000000
> M 56.666667
> Name: English, dtype: float64
<練習問題 1>
以下のデータに対して、Moneyが500以上の人を絞り込んで、レコードを表示させてください。
from pandas import Series,DataFrame
import pandas as pd
attri_data1 = {'ID':['1','2','3','4','5']
,'sex':['F','F','M','M','F']
,'Money':[1000,2000,500,300,700]
,'name':['Saito','Horie','Kondo','Kawada','Matsubara']}
attri_data_frame1 = DataFrame(attri_data1)
attri_data_frame1[attri_data_frame1['Money'] >= 500]
<練習問題 2>
上記のデータに対して、男女別（MF別）の平均Moneyを求めてください。
attri_data_frame1.groupby("sex")["Money"].mean()
> sex
> F 1233.333333
> M 400.000000
> Name: Money, dtype: float64
<練習問題 3>
上記のデータに対して、以下のデータの同じIDの人をキーとして、データをマージしてください。そして、MoneyとMathとEnglishの平均を求めてください。
attri_data2 = {'ID':['3','4','7']
,'Math':[60,30,40]
,'English':[80,20,30]}
attri_data_frame2 = DataFrame(attri_data2)
pd.merge(attri_data_frame1, attri_data_frame2)
2.1.4
<練習問題 1>
y = 5x + 3 (xは-10から10の値)のグラフを書いてみましょう。
def exam41func(x):
return 5 * x + 3
x = np.array([-10., 10.])
plt.plot(x, exam41func(x))
<練習問題 2>
先ほどのsin関数に加えて、cos関数のグラフも書いてください。
# 2行1列のグラフの1つ目
plt.subplot(2,1,1)
x = np.linspace(-10, 10,100)
plt.plot(x, np.sin(x))
plt.plot(x, np.cos(x))
# 2行1列のグラフの2つ目
plt.subplot(2,1,2)
y = np.linspace(-10, 10,100)
plt.plot(y, np.sin(2*y))
plt.plot(y, np.cos(2*y))
plt.grid(True)
<練習問題 3>
0から1の値をとる一様乱数を1000個、2組発生させて、それぞれのヒストグラムを書いてみましょう。結果はどうなっていますか。また、1000個だけではなく、100個や10000個などでも実施してみましょう。何かわかることはありますか。なお、それぞれのヒストグラムは別のグラフに表示するために、plt.subplotを利用してください。なお、一様乱数とは、ある数から数まで等確率で発生する乱数のことをいい、np.random.uniformを使ってください。
# histogram
random.seed(0)
plt.subplot(2,2,1)
exam33RandomArray1 = np.random.uniform(0, 1, 10**3)
plt.hist(exam33RandomArray1, bins=50,range=(0,1))
plt.subplot(2,2,2)
exam33RandomArray2 = np.random.uniform(0, 1, 10**3)
plt.hist(exam33RandomArray2, bins=50,range=(0,1))
plt.subplot(2,2,3)
exam33RandomArray3 = np.random.uniform(0, 1, 10**2)
plt.hist(exam33RandomArray3, bins=50,range=(0,1))
plt.subplot(2,2,4)
exam33RandomArray4 = np.random.uniform(0, 1, 10**4)
plt.hist(exam33RandomArray4, bins=50,range=(0,1))
plt.grid(True)
2.2 総合問題
2.2.1 モンテカルロ法
乱数を発生させる方法を使って、円周率を求めるプログラムを作成してみましょう。（なお、このアプローチをモンテカルロ法といいます）
（1）一様分布に従う乱数を2組発生させて、それぞれ10,000個の一様乱数を作ってみましょう。なお、一様乱数とは、ある数から数まで等確率で発生する乱数のことをいい、np.random.uniformを使います。使い方としては、0から1までの数を10個発生させる場合は、np.random.uniform(0.0, 1.0, 10)とします。
（2）x−y軸を使った中心(0,0)、半径1の円と、長さ1の正方形を考え、円の面積は$\pi$となり、正方形の面積は1となります。ここで先ほどのxとyの組み合わせの乱数10000個のうち、円の内部に入る点は何組あるでしょうか。ここで、円の内部に入るとは、x−y座標の原点から点 (x, y) のベクトルの長さを求め、それが1より小さくなる場合を判定基準とします。なお、その長さを求めるために、ユークリッドノルム($\sqrt{x^2 + y^2}$)を使い、 math.hypot(x,y)で計算できます。（さらに、余裕があれば、円の中に入ったxとyの組み合わせと外に出たxとyの組み合わせをプロットして図にしてみましょう。）
（3）半径1の1/4の円の面積と長さ1の長方形の面積の比は、$ \pi /4 : 1$となりますので、これと先ほどの結果を利用して、円周率を求めてみましょう。
# MatPlotLibは描写点数が多いと重くなるのがつらいですね。
import math
exam41RandomArrayX = np.random.uniform(0.0, 1.0, 10**4)
exam41RandomArrayY = np.random.uniform(0.0, 1.0, 10**4)
insideCounter = 0
for x, y in zip(exam41RandomArrayX[:], exam41RandomArrayY[:]):
if math.hypot(x, y) <= 1:
insideCounter += 1
plt.plot(x, y, 'o')
print(insideCounter, exam41RandomArrayX.shape[0])
# 四分円の面積 / 第一象限四角形の面積 = [1/4 * (1 * 1) * pi(四分円)] / [1 * 1] = 1/4 * pi
# -> pi = 4 * 四分円の面積 / 第一象限四角形の面積
montePi = 4 * (insideCounter / exam41RandomArrayX.shape[0])
print(montePi)
> 7885 10000
> 3.154
