#!/usr/bin/env python
# coding: utf-8
# Table of Contents
#
# # Breast Cancer Wisconsin (Diagnostic) Data Set
#
#
# ## Attribute Information:
#
# 1. ID number
# 1. Diagnosis (M = malignant:-1, B = benign:1) M:悪性，B:良性
# 1. 3-32
#
# Ten real-valued features are computed for each cell nucleus:
#
# * 半径radius (mean of distances from center to points on the perimeter)
# * テクスチャtexture (standard deviation of gray-scale values)
# * 境界の長さperimeter
# * 面積area
# * なめらかさsmoothness (local variation in radius lengths)
# * コンパクトさcompactness (perimeter^2 / area - 1.0)
# * くぼみ度合いconcavity (severity of concave portions of the contour)
# * くぼみの数concave points (number of concave portions of the contour)
# * 対称性symmetry
# * フラクタル次元fractal dimension ("coastline approximation" - 1)
#
# のそれぞれのmean, stderr, worst数値を保持している．
#
# http://people.idsia.ch/~juergen/deeplearningwinsMICCAIgrandchallenge.html
#
# ## 分類器
# 用意された訓練(training)データには，
# $\boldsymbol{A}$に上に記した特徴量が，
# $\boldsymbol{b}$に悪性(-1)か良性(1)かを
# 示す数値が入っている．
#
# 与えられた特徴ベクトル$\boldsymbol{y}$に対し，
# 細胞組織が悪性か良性かを分類する関数$C(\boldsymbol{y})$を選び出すプログラムを作成しよう．
# ## 仮説クラス
# 分類器は可能な分類器の集合(**仮説クラス**)から選ばれる．この場合，仮説クラスとは特徴ベクトルの空間$\mathbb{R}^D$から$\mathbb{R}$への線形関数$h(\cdot)$である．すると分類器は次のような関数として定義される．
#
# $$
# C(\boldsymbol{y}) =
# \left\{ \begin{array}{ccc}
# +1 & {\rm when} & h(\boldsymbol{y})\geq 0\\
# -1 & {\rm when} & h(\boldsymbol{y})<0
# \end{array} \right.
# $$
#
# 各線形関数$h:\mathbb{R}^D \rightarrow \mathbb{R}$に対して，
# 次のような$D$ベクトル$\boldsymbol{w}$が存在する．
# $$
# h(\boldsymbol{y}) = \boldsymbol{w}\cdot \boldsymbol{y}
# $$
# したがって，そのような線形関数を選ぶことは，結局$D$ベクトル$\boldsymbol{w}$を
# 選ぶことに等しい．特に，$\boldsymbol{w}$を選ぶことは，仮説クラス$h$を
# 選ぶことと等価なので，$\boldsymbol{w}$を**仮説ベクトル**と呼ぶ．
#
# 問題を単純化すると，分類器を単なるベクトルとみなして，データとの掛け算で予測がつきます．本来は予測を-1,1とかに投影しないといけないんですが，単純化のためにそのままの値を用います．問題はどうやってこの仮説ベクトル$\boldsymbol{w}$の各要素の値を決定するか？ですよね．
#
# # 最急降下法
#
# 損失関数に
# $$
# L(w)=\sum_{i=1}^n (A_i \cdot \boldsymbol{w} - b_i)^2
# $$
# を選ぶと，ベクトル$\boldsymbol{w}$の
# $j$偏微分は，
# $$
# \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}$は$A_i$の$j$番目の要素を意味する．
# この偏微分$\frac{\partial L}{\partial w_j}$を
# $\boldsymbol{w}_j$の勾配(slope)として，$L(w)$の極小値(local minimum)を求める．
#
# このような探索方法を最急降下法(steepest descent method)と呼ぶ．
# ## print_w
#
# 出てきた$\boldsymbol{w}$の$j$要素をきれいに表示する関数を用意しておきます．
# In[16]:
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[17]:
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
# In[18]:
A[0]
# In[19]:
b[0]
# ## 最急降下法によるw探索(steepest descent)
# In[20]:
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
print_w(w)
# # 結果
# In[21]:
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[22]:
show_accuracy(A, b, w)
# In[23]:
tmp = np.fromfile('./validate_A.data', np.float64, -1, " ")
A = tmp.reshape(260,30)
tmp = np.fromfile('./validate_b.data', np.float64, -1, " ")
b = tmp.reshape(260,1)
show_accuracy(A, b, w)
# # QR decomposition
#
# QR分解を使うとより簡単に最小値を求めることができる．
# 行列$A$は正方行列でないので，逆行列をもとめることができない．
# しかし，その場合でも$||A.w -b ||^2$を最小にする$w$を求めることができる．
#
# QR分解によって，$n \times m$行列は
# $$
# A = QR
# $$
# と分解される．ここで，Qは$n \times m$行列，Rは$m \times m$の正方行列で，逆行列を求めることができる．
#
# $|| Aw - b ||$が最小となるのはQRを使って，
# $$
# Q.R.w=b \\
# R.w = Q^t.b \\
# R^{-1}.R.w = R^{-1}.Q^t.b
# $$
# となりそう．
# In[24]:
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)
q, r = np.linalg.qr(A)
# In[25]:
ww = np.linalg.inv(r).dot(np.transpose(q).dot(b))
# In[26]:
q.shape
# In[27]:
print(r[0,0:5])
# In[28]:
show_accuracy(A, b, ww)
# In[29]:
print_w(ww)
# In[30]:
tmp = np.fromfile('./validate_A.data', np.float64, -1, " ")
A = tmp.reshape(260,30)
tmp = np.fromfile('./validate_b.data', np.float64, -1, " ")
b = tmp.reshape(260,1)
show_accuracy(A, b, ww)