#!/usr/bin/env python
# coding: utf-8
#
Table of Contents
#
#
#
#
# 数式処理group work-3(線形代数)解答例
#
#
#
# file:/~/python/doing_math_with_python/symbolic_math/group_works/group_work_3_ans.ipynb
#
# cc by Shigeto R. Nishitani 2009-2018
#
#
# ``` ruby
# ruby ../bin/pick_works_from_ans.rb gw3_ex_ans.ipynb -1 '' '14'
# ```
# # データ読み込み
#
# train_A.dataには特徴量が,train_b.dataには悪性(1)か良性(-1)かを示す数値が入っている.
# 訓練(tain(ing))データを読み込み,
# 仮説ベクトル$\boldsymbol{w}$の初期値を全て0.0001として,
# 最初の30データの正誤を表示せよ.
# テキストのshow_accuracyを少し改良すればできる.
#
# In[1]:
import numpy as np
tmp = np.fromfile('./train_A.data',
np.float64, -1, " ")
A = tmp.reshape(300,30)
tmp = np.fromfile('./train_b.data',
np.float64, -1, " ")
b = tmp.reshape(300,1)
w = np.zeros(30).reshape(30,1)
for i in range(30):
w[i] = 0.0001
# In[2]:
def print_w(w):
params = ["radius", "texture","perimeter","area",
"smoothness","compactness","concavity",
"concave points",
"symmetry","fractal dimension"]
print(" (params) : ",end="")
print(" (mean) (stderr) (worst)")
for i, param in enumerate(params):
print("%18s:" %param, end="")
for j in range(3):
print("%13.9f" % w[i*3+j], end="")
print()
# In[3]:
print_w(w)
# In[4]:
def show_accuracy(mA, vb, vw):
# M:悪性(-1),B:良性(1)
correct,safe_error,critical_error=0,0,0
predict = mA.dot(vw)
n = vb.size
for i in range(n):
if predict[i]*vb[i]>0:
correct += 1
elif (predict[i]<0 and vb[i]>0): # 良性なのに悪性と予測:再検査
safe_error += 1
elif (predict[i]>0 and vb[i]<0): # 悪性なのに良性と予測:見落とし
critical_error += 1
print(" correct: %4d/%4d" % (correct,n))
print(" safe error: %4d" % safe_error)
print("critical error: %4d" % critical_error)
# In[5]:
def show_first_n_data(mA, vb, vw, nn):
# M:悪性(-1),B:良性(1)
correct=0
predict = mA.dot(vw)
n = nn
for i in range(n):
print("%4d predict: %10.5f" % (i, predict[i]))
print("M(-1)or B(1): %4d" % vb[i])
if predict[i]*vb[i]>0:
correct += 1
print(" correct: %4d/%4d" % (correct,n))
# In[6]:
show_first_n_data(A, b, w, 30)
# # 距離
# 1. 行列$A$, ベクトル$w,b$の形状を確かめよ.
# 1. また,$A.w$の形状を確かめよ.
# 1. さらに$A.w-b$の距離の2乗$||A.w -b ||^2$を計算せよ.
# 1. $A.w$と$b$の距離とは乳がんの分類器においては何を意味するか?
# In[7]:
print(A.shape)
print(w.shape)
print(b.shape)
print(np.dot(A, w).shape)
LL = np.dot(A, w) - b
print(np.dot(LL.transpose(),LL))
# A[i]にはデータが入っている.またwは分類器である.これらの積$A[i].w$は悪性か良性かを判断
# する数値を与える.従ってそのベクトルの距離は,全てのデータに対する正誤の2乗和となる.
#
#
# # 係数ベクトルdLw
# 最急降下法による仮説ベクトル$\boldsymbol{w}$の最適化を試みる.
# 最急降下法の概念図を以下に示した.損失関数の値$L(\boldsymbol{w})$
# をz軸にとって(x,y)平面を$\boldsymbol{w}$と見立てて,
# その勾配$dL/dw$に従って極小値を求めるステップを
# 刻んでいく様子を示している.
# 単なるイメージ図なんで,コードの中身は無視してください.
#
# In[8]:
get_ipython().run_line_magic('matplotlib', 'notebook')
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
def f(x,y):
return x**2+y**2+x*y
x = np.arange(-np.pi, np.pi, 0.02)
y = np.arange(-np.pi, np.pi, 0.02)
X, Y = np.meshgrid(x, y)
Z1 = f(X,Y)
x_p = [2, 1, 0.5]
y_p = [2, 1, 0.5]
z_p = [f(x_p[0], y_p[0]),f(x_p[1], y_p[1]),
f(x_p[2], y_p[2])]
fig = plt.figure()
plot3d = Axes3D(fig)
plot3d.plot(x_p, y_p, z_p, "o-", color="red")
plot3d.plot_surface(X,Y,Z1,alpha=0.6)
plt.show()
# 損失関数の偏微分
# $$
# \begin{aligned}
# \frac{\partial L}{\partial w_j} &=
# \sum_{i=1}^n \frac{\partial}{\partial w_j}(A_i \cdot w -b_i)^2 \\
# &= \sum_{i=1}^n 2(A_i \cdot w -b_i) A_{ij}
# \end{aligned}
# $$
# の最後の式の$A_{ij}$の係数ベクトルをdLwとして求めよ.その次元をshapeで確かめよ?
# In[9]:
dLw = A.dot(w)-b
print(dLw.shape)
# # wの更新
#
# 係数ベクトルdLwとAのdot積が勾配ベクトルとなる.
# $$
# w = w - \sigma(dLw^t \cdot A)^t
# $$
# として仮説ベクトル$w$を勾配に従って進めたベクトルを求めよ.
# ここで$\sigma$はステップ幅と呼ばれ,勾配に従ってどの程度進むかを調整するパラメータで,大きすぎると最適値を通り越し,小さすぎると最適値にたどり着くまでに繰り返し(iteration)が多くなる.ここでは,`3.0*10**(-9)`程度とせよ.
# In[10]:
loop, sigma = 300, 3.0*10**(-9)
print_w(w)
LL = np.dot(A, w) - b
print(np.dot(LL.transpose(),LL))
w = w - (dLw.transpose().dot(A)).transpose()*sigma
print_w(w)
LL = np.dot(A, w) - b
print(np.dot(LL.transpose(),LL))
# # 最急降下の繰り返し
# 先ほどの漸近操作を300回程度繰り返し,その前後でwを表示してみよ.
# 最初の30データの予測値を比較せよ.
# In[11]:
loop, sigma = 300, 3.0*10**(-9)
for i in range(loop):
dLw = A.dot(w)-b
w = w - (dLw.transpose().dot(A)).transpose()*sigma
LL = np.dot(A, w) - b
print_w(w)
show_accuracy(A, b, w)
show_first_n_data(A, b, w, 30)
# # QR分解
#
# 行列$A$のQR分解を行い,Q, R行列の次元をshapeで確かめよ.
# In[12]:
q, r = np.linalg.qr(A)
print(q.shape)
print(r.shape)
# # 結果
# 仮説ベクトル$\boldsymbol{w}$の最適値
# $$
# ww = R^{-1}.Q^t.b
# $$
# を求めよ.その値と精度を確かめよ.また距離の2乗
# $$
# ||A.w -b ||^2
# $$
# が下がっていることを確かめよ.
# In[13]:
ww =np.linalg.inv(r).dot(np.transpose(q).dot(b))
print(ww.shape)
print_w(ww)
show_accuracy(A, b, ww)
# In[14]:
LL = np.dot(A, ww) - b
print(np.dot(LL.transpose(),LL))
# In[15]:
show_first_n_data(A, b, ww, 30)
# In[ ]:
# In[ ]: