読者です 読者をやめる 読者になる 読者になる

実体と情報のはざま

いつかデータサイエンティストになって世界を読み解く仕事がしたい!

Irisデータセット事始め

 データサイエンスの世界では有名なIrisデータセット(花の形状データ4種を変量とする150個の個体とその種類が入っている。)を使ってデータサイエンティストっぽいことをやってみたいと思い、少しデータをいじってみた。しかし、まだまだPythonとそのエコシステムが全く身についていない自分に愕然とした。。今回やったことはエクセルでやれば15分もあればできると思われる作業にすぎない。しかし、多分5時間ぐらいはかかってしまった~。
読み込みで苦戦
 今回使ったIrisのデータは最近アカウントを取得したKaggleからダウンロードしたのだが、本やネットによるとAnacondaにも同梱されているとのこと。しかし、同梱されているファイルから読み出せるのは限られたデータであって、通常は外部からCSVとしてゲットするものだ。私は汎用的な方法を身に着けたいと考えているので、今回はダウンロードしたCSVファイルを使った。
f:id:myuteru:20170528234950p:plain
 しかし、ここで困ったことが。デフォルト設定でNumpyの配列として読み込もうとしたところ、テキスト部分が読み込めない。しょうがないのでデータを文字列として読み込んだ後、数値部分はもう一度フロートに変換しなおして使うことにした。要するに、色んなデータ型が混在したNumpy配列を作れなくて苦労したって話。。いまだに分からないが時間切れということで、手持ちの技術で進める。
Pyplotには慣れてきた!
 とりあえず読み込みは完了したので、そのままプロットして様子を見る。Pyplotに関してはちょっと前に、できたのできないのでギャーギャーしながらもなんとか勝手が分かってきたので、今回はすんなりプロットできた。
まずは、生データをID順にプロットしてみた。
f:id:myuteru:20170529000542p:plain
同様に、標準化したデータをプロット。
f:id:myuteru:20170529001000p:plain
 標準化は過去にブログで書いた気がする。たった一行で標準化ができちゃうなんて、行列が直接扱えるPythonならではだな。グラフから、データセットに50個おきに段差があるということかな。概観するってデータ処理には大事なことだから、些細な情報も大切にしなくちゃね。
Pandasとseabornって!?
 次にやってみたのは、データ間の相関関係を概観するための散布図行列の作成。そのスマートなやり方が、こないだ買った本に出ていていいなと思って。仕事で似たようなことをエクセルを使ってやったことがあるんだけど、まあ、時間がかかって仕方がなかった。グラフ枚数が増えるとちょっとバグってきてたしね。それがちょいちょいってできるんなら嬉しいなと。
 だが、簡単ではなかった。まず、Pandasのデータフレームなるものの存在になかなか気づかず、seabornとの連携がうまくできずにいた。どこが難しかったかというと、軸にラベルを入れるやりかたが分からなかった。でも、最初からPandasのデータフレームでデータを読み込んでいれば簡単だった。軸のラベルが自動的についてくるから。ただ、Pandasによる読み込みがどこか邪道な気がしてしまい今回はパス。
 そんなこんなでできたのがコレ。
f:id:myuteru:20170529003122p:plain
 実はseabornのギャラリーに同じものがある。有名なデータセットだとこういうことが起きるんだね。このグラフを眺めてみると、PetalLengthとPetalWidthに相関関係がありそうなことが見て取れる。花びらの長さと幅ね。今回はお試しだからIris(アヤメ)という花の実体について詳しく知ろうとはしないけど、本来は情報が実体の何を示しているのかを注意深く考察することが大切だと思うな。
 一連のグラフは下記プログラムで生成した。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

Iris=np.genfromtxt('Iris.csv',delimiter=',',dtype=str)#まずは全部読み出し
X=np.copy(Iris[1:,1:5])#データ行列を抜き出し
X=np.asfarray(X)#型変換str→float64
Y=np.copy(Iris[1:,5])#品種の列
n=X.shape#データ行列数(行,列)
c=np.array(['b','c','r','g'])#プロット用の色を準備
#生データのプロット
fig1=plt.figure(1)
for j in range(0,n[1]):
    plt.scatter(np.arange(0,n[0]),X[:,j],
                    c=c[j],alpha=0.6,s=50,label=Iris[0,1+j])
plt.legend(loc=2)#左上に凡例を表示
plt.xlabel('Id')#X軸のラベルを指定
plt.ylabel('std data')#Y軸のラベルを指定
#標準化したデータのプロット
fig2=plt.figure(2)
for j in range(0,n[1]):
    X[:,j]=(X[:,j]-X[:,j].mean())/X[:,j].std()#データの標準化
    plt.scatter(np.arange(0,n[0]),X[:,j],
                c=c[j],alpha=0.6,s=50,label=Iris[0,1+j])
plt.legend(loc=2)#左上に凡例を表示
plt.xlabel('Id')#X軸のラベルを指定
plt.ylabel('std data')#Y軸のラベルを指定

df=pd.DataFrame(X)#pandasのデータフレームに変換
df.columns=[Iris[0,1],Iris[0,2],Iris[0,3],Iris[0,4]]#カラム名を付ける
g=sns.PairGrid(df)#散布図行列
g=g.map(plt.scatter,edgecolor='w',s=30,color='g')

 今回は良く知られているデータをよく知られた手法で概観した。次は機械学習をやってみたい。今週末は嫁と子供が実家に遊びに行ったりして勉強する時間が取れたので、ニューラルネットワークの自動微分とか勉強できたし。下記の本にはPythonを用いてディープラーニングの実装を丁寧に分かり易く書いてある。自動微分という文言は用いてないが、計算グラフを用いた誤差逆伝播法という形で解説している。また、活性化関数や損失関数を実装する際の簡潔なコードや注意点が書かれていて、初心者(私)やライブラリを使わないでやりたい人には良いと思う。ReLU関数が一行で定義できることとかソフトマックスの実装時の注意点を守らないとまともな値が出力できないこともあるなど、事前に知っておいて良かったと思える事柄を学べた。

 *アフィリエイトの画像が貼ってあるブログってカッコイイと前から思っていて申し込みをした。ついにアマゾンのリンクが貼れて嬉しい!お気に入りの本も紹介できて嬉しい!会社関係含めて身近に紹介できる人なんていないし。。
今日は、ここまで!

Kaggleにアカウント

 Kaggleのアカウントを取得した。相変わらずの英語力なので読むのは遅いが、エキスパートたちのやりとりが面白くてついつい投稿を読んでしまう。。
 サインインして初めにやったことは、タイタニックの件のチュートリアルの読み込みとデータセットのダウンロード。日本のサイトにもこれに挑戦した様子などを自分でレポートした記事が何件かあり、私がその記事を読んで大体の内容は知ってしまっている。なので、謎解きというか課題解決的な面白さは全くない。でも、実際にデータを手にしてみるとやっぱり嬉しい~。解析とかはまだやってないけど、徐々にやっていくつもり。
 もう一つダウンロードしたものが、Irisのデータセット。アヤメの花の形状データから品種を分類するという有名な課題。コンペのレベルとは程遠いけど、有名なデータセットにはすべて触れてみたいし、比較的簡単そうでいい練習になるかなと。
 金曜日を待たずして仕事とデータサイエンスの勉強疲れが出てきた。
 今日はここまで!

Kaggleっていい!

 そのうちKaggleのコンペに挑戦したいと思ってる。最近はちょっと寄り道っぽくSOMにはまってみたりしたけど、Pythonの扱いにも少しづつ慣れてきたし、そろそろKaggleのサイトに入り浸ってみようかな~と思っている。
 だが、その前に若干問題が。。それは私の英語力!Kaggleに投稿されている文章が読めない!正確には、読めるけど遅い。。楽しみながら挑戦したいと思っていて、「最短コースで上位に食い込んでやる」みたいな考えは全くないが、実力者の考え方は学びたい。仕事で時々USP(米国特許)読んだりしてるし、TOEICスコアもギリギリ海外勤務可能レベルなんだけどね。実践レベルで使える英語が身についてないってことかな。
 Kaggleのサイトでは有名なタイタニックの生き残りの話とかすっごい面白いと思うし、もしコンペで一番になれたら賞金ももらえるっていう夢があるし、挑戦し甲斐に満ち溢れている。仕事と子育てに追われながら日々を過ごしているおじさんにはハードルがべらぼうに高そうだけどね~~。不思議とやる気だけはある。
 もともと物理学が好きで仕事も物理系の仕事をしてきたんだけど、ちょっと基礎物理から外れるとすぐモデルが通用しなくなるし、工夫の余地というかクリエイティブな仕事というかエキサイティングな部分があんまりない。工業的においしい分野は物理屋的には面白くないのが普通。かといって超対称性理論とかの先端理論もかじってみたけど実体と情報のはざま(実験困難なスケールと抽象的すぎる理論のはざま)に立たされてどちらにもすすめないもどかしさを感じた。おっとブログのタイトルがでちゃったよ。。
 で、今話題で物理のお隣さんみたいな分野といえばデータサイエンスですよね。これがしっくりきた。数学モデルの扱いは得意な方だしプログラミングもかじったことあるし。こういう流れでこの分野に入っていく人多いと思うな。特に今の若い人たちは学校で普通にPC学んでくるし。
 私が初めてPCに触れたのは初代ファミコンにはまった小学生時代になんとか買ったMSXというPCだった。チャンネルが手回しの15inchそこそこのブラウン管テレビにつないで、ベーシックでスプライト(背景と独立して動く部分)を動かしてゲームっぽいやつを作ったっけ。テレビを独占すると親がうるさくて当時は存分にできなかったな~。大人になった今こそやってやるぞ!妻に怒られない範囲で。。
 今日はここまで!

入力データを近似するSOM

 前回は、自己組織化マップで色分けをやってみた。今回は、「マップ」の方ではなく入力ベクトルと参照ベクトル(プログラム中の行列Mでどんどん変わる)の「ベクトル空間」を可視化してみたいと思う。色分けは、入力ベクトルよりも参照ベクトルの方がずっと多いのだが、今回は入力ベクトルの方が多い。この場合、参照ベクトルは入力ベクトルを近似するというか寄り添うように分布することになる。
時々エラー。でも一応できた!
プログラムはこんな感じ。入力は円。

import numpy as np
import matplotlib.pyplot as plt
########################パラメータの設定
n,f=50,2#入力ベクトル数,#入力ベクトルの要素数
sgm0=0.01#誤差に関するパラメータ
T=100#繰り返し回数
a=5#時定数
cx,cy=3,2#マップのx,y方向の数
###########################行列の設定
X=np.zeros((n,f))#入力行列
for i in range(0,n):
    X[i,0]=np.cos(2*np.pi*i/(n-1))
    X[i,1]=np.sin(2*np.pi*i/(n-1))
fig1=plt.figure(1)
plt.scatter(X[:,0],X[:,1],color='b',s=20)
M=np.random.random((cx*cy,f))*0.3-0.15#マップ中の各サイトの値が入った行列
plt.scatter(M[:,0],M[:,1],color='m',s=500,alpha=0.2)
R=np.zeros((cx*cy,2))#マップ中のサイトの位置が入った行列
for y in range(0,cy):
    for x in range(0,cx):
        k=x+(cx*y)
        R[k,0]=x
        R[k,1]=y
Ita=np.zeros((n,cx*cy))#更新に使う係数
###########################繰り返しの実行
for t in range(0,T):
    sgm=sgm0*np.exp(-t/a)
    for i in range(0,n):
        min=1
        for j in range(0,cx*cy):
            d=np.linalg.norm(X[i,:]-M[j,:])
            if d<min:
                min=d
                cj=j
        for j in range(0,cx*cy):
            D=np.linalg.norm(R[j,:]-R[cj,:])
            Ita[i,j]=np.exp(-((D/sgm)**2)/2)
    s=np.sum(Ita,axis=0)#列方向に和をとる(axis=1なら行方向)
    for j in range(0,cx*cy):
        p=np.zeros((f))
        for i in range(0,n):
            p=p+(Ita[i,j]*X[i,:])
        M[j,:]=p/s[j]
plt.scatter(M[:,0],M[:,1],color='r',s=500,alpha=0.5 )

 結果がこれ。うまくいったケースのもの。入力ベクトルが形作る円を少数の参照ベクトルが近似している様子が見られる。円の中心付近のマゼンタ色で散らばっている点が初期の参照ベクトルで、繰り返し計算の末が赤い点である。プログラムの未解決課題は時々「ゼロで割る」ことによるエラーが出てしまうこと。。前回色分けの時に収束させる関数を変えたんだけどあんまり関係なかったみたい。
f:id:myuteru:20170521000210p:plain
 次回は何かしらを”分類”してみたいな~。SOMからいったん離れて”分類”手法の古典から学びたいな。
今日はここまで!

自己組織化マップに注力

 前回は、初めて自己組織化マップ(SOM)をPythonで書いてみたが、いざやってみると自分が根本をいまいち理解できていないことがはっきりした。今も実はまだちゃんと理解できていない。ネットや本で調べてみたが様々な考え方があるようだし、根本的な部分で腑に落ちるものはなかった。ただ、平均寿命の半分を生きてきた経験から、「分からない状態をキープする」ことで自分が成長できるのも知っている。SOMについてはこの状態にしておきたい。というわけで、もうしばらくはSOMの周辺に住み込んでいたいと思う。
 理屈はともかくSOMらしいものはできた!
 前回のプログラムはいまいち納得できないことばかりだったので書き直した。一番重要な変更点はこの部分のσの関数形。
f:id:myuteru:20170519231220p:plain
これは、マップ上において着目地点の値の更新の影響をどれくらい遠くまで与えるかを表現した式。今回、関数系を変えたのは、あとで出てくる計算のゼロで割ることのエラー回避のため。減少関数ならいろんな選択肢があるが、シンプルで挙動が分かり易いのを採用した。あとは、プログラム構造のダサい感じを少しでもスマートになるように工夫した。(けど、大差なしか~。。)
プログラムはこうなった。

import numpy as np
import matplotlib.pyplot as plt
########################パラメータの設定
n,f=5,3#入力ベクトル数,#入力ベクトルの要素数
sgm0=10#誤差に関するパラメータ
T=50#繰り返し回数兼時定数
cx,cy=20,15#マップのx,y方向の数

###########################行列の設定
X=np.zeros((n,f))#入力ベクトル行列
X[0,:]=([1,0,0])#R
X[1,:]=([0,1,0])#G
X[2,:]=([0,0,1])#B
X[3,:]=([0,1,1])#?
X[4,:]=([1,1,0])#?
M=np.random.random((cx*cy,f))#マップ中の各サイトの値が入った行列
R=np.zeros((cx*cy,2))#マップ中のサイトの位置が入った行列
for y in range(0,cy):
    for x in range(0,cx):
        k=x+(cx*y)
        R[k,0]=x
        R[k,1]=y
Ita=np.zeros((n,cx*cy))#更新に使う係数
fig1=plt.figure(1)
for j in range(0,cx*cy):
    plt.scatter(R[j,0],R[j,1],color=(M[j,0],M[j,1],M[j,2]),s=100,marker='s')

###########################繰り返しの実行
for t in range(0,T):
    sgm=sgm0*np.exp(-t/T)#
    for i in range(0,n):
        min=1
        for j in range(0,cx*cy):
            d=np.linalg.norm(X[i,:]-M[j,:])
            if d<min:
                min=d
                cj=j
        for j in range(0,cx*cy):
            D=np.linalg.norm(R[j,:]-R[cj,:])
            Ita[i,j]=np.exp(-((D/sgm)**2)/2)
    s=np.sum(Ita,axis=0)#列方向に和をとる(axis=1なら行方向)
    for j in range(0,cx*cy):
        p=np.zeros((f))
        for i in range(0,n):
            p=p+(Ita[i,j]*X[i,:])
        M[j,:]=p/s[j]

fig2=plt.figure(2)
for j in range(0,cx*cy):
    plt.scatter(R[j,0],R[j,1],color=(M[j,0],M[j,1],M[j,2]),s=100,marker='s')

 結果は、こんな感じ。
初期はランダム。
f:id:myuteru:20170519232624p:plain
繰り返し計算で自己組織化した状態がこれ。グラデーションがきれい。
f:id:myuteru:20170519232644p:plain
入力として5色をRGBで指定した。
 プログラムがうまく動いたので調子に乗ってマップのサイズ(プログラム中の記号ではcx,cy)を大きくしてみると、急に計算時間が長くなった。ビッグデータをSOMで分類している方々は何かしらの高速化の仕組みを持っているのだろうな。
 次回は、もっと”分類”や”次元圧縮”といった側面に触れたい~。
今日はここまで!

自己組織化マップSOM(Self-OrganizingMap)に挑戦した!

 今回は、ずっと気になっていたSOMに挑戦!アルゴリズムはコーエン博士の本を購入してなんとか分かってきていたので、さっそく勉強中のPythonでプログラミングを試してみた。ライブラリーを使えば簡単なのだろうけど、いつか自分でこの手のアルゴリズムを作りたいと思っている(何年後?)ので、フルスクラッチでやってみたい。
変な結果に。。
 3日がかり。。といっても仕事と子育てで実働6時間弱ってとこだけど。
 やってみたのは、自己組織化マップの世界ではおなじみ(?)の色の分類。はじめランダムなんだけど、徐々に入力した色を中心とした分類が見えてくるってやつ。何に苦労したかというと、そもそもの理論とPython。全部か。いざプログラミングを始めてみたら、理解したと思っていたバッチ方式のSOMの理論に自信がなくなってきて上記の本に何度かもどった。Pythonの行列の扱いもまだ体に染み付いてない。。行列の列の和をもとめてベクトル化する方法とか知らんかったし。
 で、なんとか書き上げたのがコレ!for文とか多いしdef(定義)も使ってなくて見通しが悪いし、素人の私が見てもブサイクなコードだとすぐわかる。。次回はもっとブラッシュアップしてキレイに仕上げたいな~。

import numpy as np
import matplotlib.pyplot as plt

########################パラメータの設定
n=3#入力ベクトルの数
f=3#入力ベクトルの要素数
sgm0=2#誤差に関するパラメータ
T=10#繰り返し回数
cx=12#マップのx方向の数
cy=12#マップのy方向の数

###########################行列の設定
X=np.zeros((n,f))#入力行列
X[0,:]=([1,0,0])#R
X[1,:]=([0,1,0])#G
X[2,:]=([0,0,1])#B
M=np.random.random((cx*cy,f))#マップ中の各サイトの値が入った行列
R=np.zeros((cx*cy,2))#マップ中のサイトの位置が入った行列
j,k=0,0
for i in range(0,cx*cy):
    j=j+1
    R[i,0]=j-1
    R[i,1]=k
    if j > cx-1 :
        j=0
        k=k+1
Ita=np.zeros((n,cx*cy))#更新に使う係数
fig1=plt.figure(1)
for j in range(0,cx*cy):
    plt.scatter(R[j,0],R[j,1],color=(M[j,0],M[j,1],M[j,2]),s=200,alpha=0.9,marker='s')

###########################繰り返しの実行
for t in range(0,T):
    for i in range(0,n):
        min=1
        for j in range(0,cx*cy):
            d=np.linalg.norm(X[i,:]-M[j,:])
            if d<min:
                min=d
                Ci=j
        sgm=sgm0*(T-t)/T
        for j in range(0,cx*cy):
            D=np.linalg.norm(R[j,:]-R[Ci,:])
            Ita[i,j]=np.exp(-((D/sgm)**2)/2)
    s=np.sum(Ita,axis=0)#列方向に和をとる(axis=1なら行方向)
    for j in range(cx*cy):
        p=np.zeros((n))
        for i in range(0,n):
            p=p+(Ita[i,j]*X[i,:])
        M[j,:]=p/s[j]
fig2=plt.figure(2)
for j in range(0,cx*cy):
    plt.scatter(R[j,0],R[j,1],color=(M[j,0],M[j,1],M[j,2]),s=200,alpha=0.9,marker='s')

この結果はこんな感じ。まずランダム。入力はRed,Green,Blueの3色。
f:id:myuteru:20170518001851p:plain
次が、自己組織化の結果。黒いところは多分エラーが出たとこ。
f:id:myuteru:20170518001911p:plain
 この結果、何かが違う。入力した色しかないし。マップを大きくするとすぐエラーが出るし。パラメーターをちょっといじるとやっぱりすぐエラーになるしで微妙~。だけど正解に近づいた気がして嬉しー!!次回はうまくできているかの検証も含めてやってみたい。
 あと、最近このブログを読んでくださる方が増えてきたようなので、プログラムの載せ方や注釈なども時間をかけたいと思う。☆マークをつけてくださった方もいらっしゃって、感無量です。今は記事を書くこと(自分のレベルアップ)に夢中で、まわりが見えてないです。はてなブログの仕組みもいまいち理解できてないです、落ち着いたら皆さんのブログを訪問したいと思います。

PythonのNumPyについて基本中の基本を学んだ。。

 ある程度Pythonにも慣れてきたし、データサイエンスっぽい何かをやってみたいと思い立った。
 まずは、まえからずっとやりたかったBL-SOM(バッチ方式の自己組織化マップ)をコーディングしようとSpyderに向かったが…書けない!配列の扱いが全然わからない!C#で画像処理をした時は、画素を行列と見立ててx[i,j]という具合に記述してforで処理を回していったが…Pythonだとどうすりゃいいのかわからない。。
というわけで、基本中の基本を勉強しました。やり方は色々あるようだけど、やっぱNumPyが一番よさそう。
 以下、忘備録。

import numpy as np

a=np.array([2,-4,1])
b=np.array([3,1,-2])
x=np.dot(a,b)#ベクトルの内積
print(x)

x=np.linalg.norm(a-b)#ノルム
print(x)

x=a-b#ベクトルの引き算
print(x)

x=np.zeros((5,3))#ゼロでできた5行3列の行列
print(x)

x=np.random.random((5,3))#乱数でできた5行列の行列
print(x)
print(x[1,1])#行列の要素(1,1)を取り出し。(角が())

x[0,0]=0#行列の要素(0,0)を0に変更
print(x)

以下、結果。
0
5.9160797831
[-1 -5 3]
[[ 0. 0. 0.]
[ 0. 0. 0.]
[ 0. 0. 0.]
[ 0. 0. 0.]
[ 0. 0. 0.]]
[[ 0.66526329 0.12998445 0.89432272]
[ 0.58101282 0.44986365 0.81146389]
[ 0.96562799 0.8162834 0.08969892]
[ 0.07098462 0.34161584 0.96946225]
[ 0.19910363 0.58869816 0.30484134]]
0.449863645793
[[ 0. 0.12998445 0.89432272]
[ 0.58101282 0.44986365 0.81146389]
[ 0.96562799 0.8162834 0.08969892]
[ 0.07098462 0.34161584 0.96946225]
[ 0.19910363 0.58869816 0.30484134]]

次回はBL-SOMに挑戦したい!また基本的なところでつまづかないか心配だが。。
今日は、ここまで!