ニューラルネットで雑な手書き文字認識

以前にもニューラルネットについて言及しましたが，かなり原初的な内容ばかりだったので今回は簡単な応用について．

C言語でニューラルネット（関数近似）

C言語でニューラルネット（2次元平面での分離）

といっても，修士の授業で作った内容そのままなので，いざ実装してみると動かないかも．

ニューラルネットの応用先として，真っ先に上がるのがパターン認識の分野でしょう．

とくに文字認識の分野ではすでにかなりの数，実用化されていると聞きます．

ここでは，最も単純な文字の認識について述べていきたいと思います．

具体的には以下の2点．

• ニューラルネット（誤差逆伝播法）の基本的な部分
• 特徴量（ニューラルネットへの入力）

パーティクルフィルタによるグリッドマップマッチング

今回，著作権的にグレーな感じがするので問題があればご一報ください．

専門ではないですが，ロボットの自己位置推定に興味を持ったのでMATLAB/simulinkでシミュレーションしてみました．
ところどころ用語が間違っていると思いますがご容赦ください．

手を付けたきっかけは，“グリッドマップのマッチングに基づく未知障害物にロバストな自己位置推定”という論文の発表を聞いたことです．

面白そうだったので今回は，これを実装したいと思います．

また，今回の記事の大部分はMyEnigmaさんの記事を参考にしています．

パーティクルフィルタやグリッドマップそれぞれについては，こちらがとてもに詳しいのですが，これらを組み合わせた実装について記載がなかったため，本記事で述べていきます．

自己位置推定は，ロボットが現在，どこに居るのかを知るための手法です．
見晴らしの良い屋外であればGPSで事足りることが多いかと思いますが，そうでない場合にはロボット自身が現在の位置を知ることは難しくなります．

自己位置推定はよく知るものだとカーナビがありますが，これは普段はGPSで位置を測定していて，トンネル内などのGPSの電波が入らないところでは，タイヤの回転数（回転角）や加速度センサの2回積分（$\simeq$位置）を使って，自己位置を推定しています（多分）．

カーナビでは，トンネル内であまり自己位置の精度が要求されない（トンネル内での分岐は少ない）ことと，道路自体がカーナビ内部にあるので推定された位置が多少ズレていても道路に引き寄せられること，最終的な運転は人間が行うことから，あまり自己位置推定は難しくない/重要でないのでは，と以前は思っていました．

ロボットの場合だと，精度の高い自己位置の推定が要求されるため，いろいろな手法が開発されてきたとのことです．あまり詳しくないのでこの辺りは省略します．

今回はその中でもパーティクルフィルタによるグリッドマップマッチングについてシミュレーションをしてみました．

C言語でニューラルネット（2次元平面での分離）

前回に引き続き，ニューラルネットワークです．

前回は，ニューラルネットによる関数近似を行いましたが，今回はちょっとパターン認識寄り．
というのも，2次元平面上のある点が，AとBどちらに分類されるのかという問題なためです．

たとえば，以下の様な状況の場合．

これを2種類に分類したいとします．
理想的な場合，このように分類できると想像できると思います．

先ほどの4点を与えたとき，大体の人はこのように線を引いて分けると思います．
このような1本の線では分離ができない問題のことを，線形分離不可能な問題と言ったりします．

線形分離可能な例はこんな感じ．

1本の直線で，分類できています．

線形分離可能な問題は，単純パーセプトロンなどでも解けますが，線形分離不可能の場合には，隠れ層を持つようなニューラルネットワークの構成である必要があります．

このような隠れ層を持つニューラルネットワークの学習方法のひとつが，誤差逆伝播法（Back Propagation）です．

というわけで，2次元での分類を行っていきます．

まず，2次元の入力（x軸，y軸）に対して出力（●，×）があるので，入力次元は2，出力次元は1です．

data.dat		rho.dat
0	0	1
0	1	0
1	0	0
1	1	0

まずは，線形分離可能なパターンから試してみます．
ここで，data.datの左側がx軸，右側がy軸で，rho.datの0が●，1が×を表しています．
よって，このデータは先ほどの直線で分離できるデータを数字で表したものとなっています．

これをニューラルネットワークに入力として与えると，以下のようになります．

このように，確りと分離できていることがわかります．
プログラム中では，positive.datに0.5以上（厳密には超過）の出力群を，negative.datに0.5以下の出力群をファイル出力しています．

次に，線形分離不可能の場合．

data.dat		rho.dat
0	0	1
0	1	0
1	0	0
1	1	1

これは最初の図に対応します．
このときの出力は，以下のようになります．

左下と右上で，×になっています．
ここからわかるように，非線形の分離問題（EX-OR問題）に対しても分離を行えていることがわかります．

プログラムリストは，前回の物から少し変えてあります．（主にTest関数内）
以下リスト．

C言語でニューラルネット（関数近似）

C言語で関数近似を行うためのニューラルネットワーク（誤差逆伝播法，B ack P ropagation N eural N etwork）のプログラムです．

プログラム中で，Input_sizeは学習用データの入力の長さ，Input_dimは学習用データの次元，test_sizeは検証用データの入力の長さです．

たとえば，以下のような点の集合を関数で近似したい場合を考えます．

このとき，学習用データは次のように与えられます．

data.dat	rho.dat
0	1
1	3
2	2
3	4

というわけで，上のように実行ファイルの場所に，data.datとrho.datを作成します．
見て分かる通り，data.datがx軸，rho.datがy軸を表しています（一次元の場合）．

この2種類のファイルから，Input_sizeとInput_dimを設定します．
この場合，Input_sizeは，データが横軸に4つなので，Input_size=4とします．
また，入力は1次元なのでInput_dim=1です．

また，test_sizeは，学習が終わったあとのデータ点の数で，x軸の0~3を何分割するかを決める滑らかさの値です．
理想的な場合学習が起こった場合では，test_size=4とすると，上の図の点と同じ位置に点が打たれるはずです．
ここでは，どのように近似されたかを詳細に見るために，test_size=100とします．

これらのデータから実行した結果が以下の図．（学習回数1,000,000回，test_size=100）

このように，学習データとして与えた点を全て通るような曲線が描かれています．
これが，ニューラルネットワークによる関数近似です．

先ほど与えた点の値と，現在のニューラルネットワークの出力の誤差がどんどん小さくなっていくことで，与えた点を通る曲線が描かれています．
この曲線を拡大してみると，100個の点からなる直線の集合です（=test_sizeの値）．

test_sizeを小さくするとこんな感じに．（学習回数1,000,000回，test_size=10）

また，学習回数が少ない場合など．（学習回数10,000回，test_size=100）

学習回数が少ないと，与えた点との近似率が落ちています．（点を通っていない）

学習回数が更に少ない場合．（学習回数1,000回，test_size=100）

ここまでくると，もはや点を通ってすらいないです．

このように，ニューラルネットワークでの関数近似は，与えられた点と現在のニューラルネットワークによる出力の誤差が小さくなるように繰り返し学習することで，与えられた点を通る出力を得られるというものです．

というわけで以下リスト．

プログラム内の変数の詳細についてまとめておくと以下のようになっています．

label	details
Input_size	入力の長さ（大きさ）
Input_dim	入力の次元
Hidden_dim	中間層の数
alpha	シグモイド関数の傾き（後述）
eta	学習係数（他サイト参照）
gain	ゲイン（後述）
test_size	学習後の検証用データの量（大きさ）
th	未使用

ここで，ニューラルネットワークに使用しているシグモイド関数は，0~1の範囲の実数しか取ることができないので，このままでは0~1の出力の範囲の点の近似しかできません．

このため，シグモイド関数にgainを掛けることで，擬似的にシグモイド関数の出力範囲を広げています．
たとえば，gain=100の場合では，シグモイド関数の出力範囲は-50~50の実数となります．
一応，負の値も取れるように，負方向にgainの半分だけずらしています．

これらの処理により，ニューラルネットワークで-50~50までの範囲の点を近似することができるようになっています．

Euler法による一階常微分方程式の解法（CR回路を例に）

少し前の予告通りEuler法による微分方程式の求解を。

まず、微分方程式というのは方程式内に微分の形を含むもので、
例えば速度のみが与えられる場合にも微分方程式となると思います。
（勝手な定義してます。ごめんなさい。）

この場合、

であり、両辺を積分することで

が得られます。

これが一般的な不定積分ですけど、数値解析の場合にはdxとdtを
「無限に小さな刻み幅」から少しゆるめて、
「有限の小さな値」として取り扱います。

有限になるということは普通の値と同様に扱うことができるようになります。

よって上の式を

とできます。そしてこれは

と変形できます。

ここでvは定数、Δtは任意の刻み幅なので微小時間Δtにおける変位、
x(t+Δt)-x(t)が求められます。x(t)は現在の位置です。
最後の式は微小時間が経過したあとの位置です。
一瞬前の場所x(t)に微小時間と速度を掛けたものを足した値が新しい位置となっています。


これを基本にして次の二階微分方程式を解きます。
画像多めにつき続きからどうぞ。

Newton-Roaphson法による方程式の求解

関数のグラフ化により数値解の個数とおおよその解が判別出来たところで、
その解を求めていく作業に入ります。 

全部書くのは大変なのでwikipediaさんのお力を借ります。
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

ニュートン法はwikipediaにもあるように、 まず、任意のx座標（wikiではxn）を決めて、
そこのy値での接線を引きます。



接線を引くと接線の傾きが0でなければ、x軸のどこかと交差します。（wikiではxn+1の点） 
そのxn+1の点を再びxnとして用い、そこのy値の接線を引きます。
そしてまた接線との交点をxn+1とします。 これの繰り返しで方程式の近似解が得られます。
x軸と接線との角をθとすると、tanθは

   

ここで右の二式で求まっていないものはxn+1だけなのでそれを変形すると、

 

となるので現在分かっているxn,f(xn),f'(xn)から未知の点xn+1を求められます。
 これを繰り返すと設定した点xnから大体一番近い解に収束していきます。

 Newton-Roaphson法を使うにあたって必要なこと。

 1.方程式が判っている。（関数のグラフ化では0.5X^3 + 0.75X^2 - 2X - 2）
 2.解の個数が判っている。（同じく3個）
 3.接線の方程式（方程式の一階微分）が判っている。（同じく1.5X^2 + 1.5X - 2）
 4.解の大体の値が判っている。

 こんな感じです。
使用時はValue of X にそこそこ解に近い適当な値（xn）を入れてください。
forの中で10回ループしてますけど、大体そのくらいで収束します。

二つ以上の解の場合、たとえば上の図だと初期値が3以下であれば解は1くらいの値に、
3以上であれば5くらいの値に収束します。
3の場合はエラーだと思います。（接線の傾きが0のためx軸と交わらない）


 以下リスト。

数学関数のグラフ表示

今回は数学関数をX-WindowSystemでグラフ化する方法を載せたいと思います。
グラフ化することによって視覚的にわかり易くなり、実数解が存在するか程度ならすぐに分かります。

ではひとまず定義。



これが表示用の関数です。

ウィンドウの大きさは480*480で、うち上下左右40pxは空白部とします。
また、数値解析における分割数は100分割で、400px内を100分割して直線で結ぶ感じです。
つまり4pxは直線だけど大きい目で見たら曲線に見えるよね、って感じかな。

線を引くためには始点と終点が必要になりますが、終点は現在の計算点、始点は前回の計算点とします。
0→1、1→2、2→3点という風になります。

以下プログラムリスト。
Linuxの人はindentコマンドで整形するなりなんなりしてくだされ。

ちなみにindentコマンドは
sudo apt-get install indent
より
indent FILENAME -OPTION
で使えます。
OPTIONは-krをつかうと幸せ。

こんな感じになります。



以下リスト。

Windowsプログラミングによるオセロ

学校の課題の休憩中にWinAPIの練習がてらいつものオセロを作ってみました。

Drawing関数はInvalidateRect()とかでWM_PAINT内で処理させた方がよかったかな。

最初に何よりも迷ったのがVisualC++でのWinプロジェクトの作成。
コンソールだとウィンドウとか出せないので、Winプロジェクトを作る必要があるんだけれど、
わざわざそんな初歩的なこと書いてないし、
なにより参考HPが基本的に以前のVerのVisualStudioだったりとか。

新規プロジェクト→Win32→Win32プロジェクトから作れます。
コンソールアプリケーションじゃだめ。

中身自体は空のプロジェクトでOK。そのあといつものように.cppを追加しとけば大丈夫。
とりあえずAPIの基本的な動き方はわかった（つもり）ので、
これからチェスとか作れればいいなぁなんて。

ちょっと前にもあったけどネタ。
http://www.gizmodo.jp/2011/06/667.html

これに対応する予定は今のところありません（笑）


以下リスト。

#include<windows.h>
#include <tchar.h>
#include<stdlib.h>
#include<stdio.h>


//const
#define WINDOW_WIDTH (430)
#define WINDOW_HEIGHT (340)


//grobal
RECT g_windowPos;
int board[10][10] = {0};
int buf_board[10][10] = {0};
int rev_board[10][10] = {0};
static int act = 1,wait = 2;//1 = player1,2 = player2
int rev_act = 2,rev_wait = 1;


//Prototype
HWND Create(HINSTANCE hInst);
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp);
void OthelloInit(HWND hWnd);
void Drawing(HWND hWnd);
int Algo(int x,int y,HWND hWnd);
void Check(int x,int y,int i,int j,int *flag);
void REV(int flag,HWND hWnd);
int NUMSTONE(int color);


//Start
int WINAPI WinMain(HINSTANCE hInst,
HINSTANCE hPrevInst,
LPSTR lpCmdLine,
int showCmd)
{
HWND hWnd;
MSG msg;


//Create Window
hWnd = Create(hInst);


if ( hWnd == NULL){
MessageBox(NULL,_T("Creating Window is failed"),_T("ERROR"),MB_OK);
return 1;
}


//Show Window
ShowWindow( hWnd, SW_SHOW);
UpdateWindow(hWnd);


OthelloInit(hWnd);


//Message Loop
while( 1 ){
BOOL ret = GetMessage( &msg,NULL,0,0);//Get Message
if( ret == 0 || ret == -1){
//if get Message of Shutdown,
//or failed GetMessage() (return -1),EndLoop
break;
} else {
//Message
TranslateMessage( &msg );
DispatchMessage( &msg );
}
}
return 0;
}


HWND Create(HINSTANCE hInst){
WNDCLASSEX wc;


// Window Class
wc.cbSize = sizeof(wc); //const size
wc.style = CS_HREDRAW | CS_VREDRAW; //style
wc.lpfnWndProc = WndProc; //window procedure
wc.cbClsExtra = 0; //Extra Info1
wc.cbWndExtra = 0; //Extra Info2
wc.hInstance = hInst; //Instance handle
wc.hIcon = (HICON)LoadImage( //Icon
NULL,MAKEINTRESOURCE(IDI_APPLICATION),IMAGE_ICON,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hIconSm = wc.hIcon; //child Icon
wc.hCursor = (HCURSOR)LoadImage( //Mouse Cursor
NULL,MAKEINTRESOURCE(IDC_ARROW),IMAGE_CURSOR,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);//Window Background
wc.lpszMenuName = NULL; //Menu Name
wc.lpszClassName = _T("Default Class Name");//Window Class Name


//Submit Window Class
if ( RegisterClassEx( &wc ) == 0){
return NULL;
}


//Window's location
g_windowPos.left = (GetSystemMetrics( SM_CXSCREEN ) - WINDOW_WIDTH )/2;
g_windowPos.top = (GetSystemMetrics( SM_CYSCREEN ) - WINDOW_HEIGHT )/2;
g_windowPos.right = g_windowPos.left + WINDOW_WIDTH;
g_windowPos.bottom = g_windowPos.top + WINDOW_HEIGHT;


//Create Window
return CreateWindow(
wc.lpszClassName, //Window Class Name
_T("Sample Program"), //Title bar
WS_OVERLAPPEDWINDOW & ~WS_THICKFRAME, //Window type
g_windowPos.left, //X
g_windowPos.top, //Y
WINDOW_WIDTH, //Window Width
WINDOW_HEIGHT, //Window Height
NULL, //Parent Window handle
NULL, //Menu Handle
hInst, //Instance Handle
NULL //Extra Data
);
}


//Window Procedure
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp)
{
static int x,y;
static BOOL draw = FALSE;
HWND hButtonWnd;
HINSTANCE hInst;
hInst = (HINSTANCE)GetWindowLong(hWnd,GWL_HINSTANCE);


switch( msg ){
case WM_LBUTTONUP:
x = LOWORD(lp);
y = HIWORD(lp);
if(x >= 20 && x <= 55+34*7 && y >= 20 && y <= 55+34*7){
Algo(x,y,hWnd);
Drawing(hWnd);
}
return 0;
case WM_RBUTTONUP:
if(act == 1){
act = 2;
wait = 1;
}else if(act == 2){
act = 1;
wait = 2;
}
Drawing(hWnd);
return 0;
case WM_CREATE:
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("UNDO"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,20,100,30,hWnd,(HMENU)1000,hInst,NULL);
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("RESTART"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,260,100,30,hWnd,(HMENU)2000,hInst,NULL);
break;
case WM_COMMAND:
switch(LOWORD(wp)){
case 1000:
REV(2,hWnd);
Drawing(hWnd);
break;
case 2000:
OthelloInit(hWnd);
return 0;
}
break;
case WM_CLOSE:
if( MessageBox(hWnd,_T("shutdown?"),_T("checking"),MB_YESNO ) == IDNO){
return 0;
}
break;
case WM_DESTROY://when Window Destroy
PostQuitMessage(0);
return 0;
}


//other Messages is DEFAULT
return DefWindowProc(hWnd,msg,wp,lp);


}


void OthelloInit(HWND hWnd){
int i,j;


act = 1;
wait = 2;
for (i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = 0;
}
}


for (i = 0;i <; 10;i++){
board[i][0] = 8;
board[i][9] = 8;
board[0][i] = 8;
board[9][i] = 8;
rev_board[i][0] = 8;
rev_board[i][9] = 8;
rev_board[0][i] = 8;
rev_board[9][i] = 8;


}
board[4][5] = 1;
board[5][4] = 1;
board[4][4] = 2;
board[5][5] = 2;
rev_board[4][5] = 1;
rev_board[5][4] = 1;
rev_board[4][4] = 2;
rev_board[5][5] = 2;


EnableWindow(GetDlgItem(hWnd,1000),FALSE);


Drawing(hWnd);
}


void Drawing(HWND hWnd){
int i,j;
char word[50];
RECT rc;
HDC hDC;
HBRUSH hBrush;
HPEN hPen;
WCHAR prev_word[50];


hDC = GetDC(hWnd);
/*
GetClientRect(hWnd,&rc);
InvalidateRect(hWnd,&rc,FALSE);
*/
hPen = CreatePen(PS_NULL,0,0);
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hPen);

SelectObject(hDC,hBrush);
Rectangle(hDC,310,160,400,200);
DeleteObject(hPen);
DeleteObject(hBrush);

//GetStockObject(BLACK_PEN);
SelectObject(hDC,(HPEN)GetStockObject(BLACK_PEN));
for (i = 0;i < 8;i++){
for (j = 0;j < 8;j++){
hBrush = CreateSolidBrush(RGB(0,255,0));
SelectObject(hDC,hBrush);
Rectangle(hDC,20+34*i,20+34*j,55+34*i,55+34*j);
DeleteObject(hBrush);
switch(board[i+1][j+1]){
case 1:
hBrush = CreateSolidBrush(RGB(0,0,0));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
case 2:
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
}
}
}
SetRect(&rc,310,100,400,150);
switch(act){
case 1:
DrawText(hDC,_T("Black turn"),-1,&rc,0);
break;
case 2:
DrawText(hDC,_T("White turn"),-1,&rc,0);
break;
}

sprintf_s(word,30,"black = %d\nwhite = %d",NUMSTONE(1),NUMSTONE(2));
mbstowcs_s(0,prev_word,30,word,_TRUNCATE);
SetRect(&rc,310,160,400,200);
DrawText(hDC,prev_word,-1,&rc,DT_WORDBREAK);
ReleaseDC(hWnd,hDC);
}


int Algo(int x,int y,HWND hWnd){
int locx,locy,i,j,flag;
for(i = 0;i < 8;i++){
if (x >= 20+34*i && x <= 55+34*i){
locx = i+1;
}
if (y >= 20+34*i && y <= 55+34*i){
locy = i+1;
}
}


if(board[locx][locy] != 0){
return -1;
}


if(locx >= 1 && locx <= 8 && locy >= 1 && locy <= 8){


/*置き石判定*/
flag = 0;
REV(flag,hWnd);
for (i = -1;i <= 1;i++){
for(j = -1;j <= 1;j++){
Check(locx,locy,i,j,&flag);
}
}
if(flag != 0){
REV(flag,hWnd);
board[locx][locy] = act;
if(act == 2){
rev_act = act;
rev_wait = wait;
wait = act;
act = 1;
} else {
rev_act = act;
rev_wait = wait;
wait = act;
act = 2;
}
}
}
return 0;
}


void Check(int x,int y,int i,int j,int *flag){
int k,num;
if(board[x+i][y+j] == wait){
num = 0;
while(1){
if(board[x+i][y+j] == wait){
num++;
} else if (board[x+i][y+j] == act){
for(k = 0;k < num;k++){
board[x-i*k][y-j*k] = act;
}
*flag = 1;
break;
} else if(board[x+i][y+j] == 8 || board[x+i][y+j] == 0){
break;
}
x += i;
y += j;
}
}
}


void REV(int flag,HWND hWnd){
int i,j;
if(flag == 0){
for(i = 0;i < 10;i++){
for(j = 0;j &t; 10;j++){
buf_board[i][j] = board[i][j];
}
}
} else if(flag == 1){
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
rev_board[i][j] = buf_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),TRUE);
} else if(flag == 2){
act = rev_act;
wait = rev_wait;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = rev_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),FALSE);
}
}

int NUMSTONE(int color){
int i,j,num;
num = 0;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
if(board[i][j] == color){
num++;
}
}
}
return num;
}