k-nnを使った簡単なパターン認識

勉強がてらに簡単にやってみました。

k-nn法(k-nearest neighbor method)は、本来確率密度を推定する際に、ノンパラメトリック（推定がパラメータに依らない）な手法をして用いられますが、これをパターン認識に応用することができます。詳しくはPRML上巻のノンパラメトリック法をご覧ください。

パターン認識と機械学習 上 - ベイズ理論による統計的予測

posted with amazlet at 11.09.30

C. M. ビショップ
シュプリンガー・ジャパン株式会社
売り上げランキング: 58427

Amazon.co.jp で詳細を見る

詳しいことは省略する（正直うまく説明できる自信がない）方向で行きます。

以下のように、各人間が支援している政党がある状況を考えます。ここで、お互いが近ければ近いほどより類似度が高い(仲が良い)とします。

ここで、真ん中にいる人が自分が支援する政党を、自分のお友達が支援している数が最も多い政党にするとします。

この時、もし仲良し3人の支援している政党を見るとすれば、民主が2人、自民が1人なので、この人は民主を支援することになります。

もし5人だったら

となります。

この考えを使って、前回の記事でも扱った手書き文字の分類をしてみます。
今回は実験を簡単にするため、文字を「０」と「１」の2種類に限定します。

流れだけ説明すると・・・

1.次元数の特徴ベクトルからなるノードから、何らかの手法を用いて類似度を計算する

2．類似度に基づいてk-nnをかける

3．正しく認識ができているか確認

・データ
16×16ピクセルで表現された手書き文字データ。各ピクセルは256階調のグレースケールで表現されている。
今回利用したデータの文字数は2115で、「１」と「０」がほぼ同数になっています。

本来は交差検定みたいなことをしないといけないと思いますが、今回は、初めからすべてのデータのk-nnを計算しておき、そのうえで各文字の近傍を見て判定をします。

まず、元のデータから各任意の文字2つの類似度を計算します。文字の集合をNとすると、N個の文字から2つの文字の類似度をすべて計算するためには、N×(N-1)/2通りの計算を行う必要があります。

プログラムは以下のような形

入力データは前回の記事で説明したような圧縮形式になっており、それを類似度の高い順にソートしたものを出力とします。類似度はコサイン類似度を利用しました。

//Dm*Dnの特徴ベクトルからcosin類似度を用いて類似度行列を作成し、降順ソートします
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <vector>
#include <algorithm>
#include<sys/times.h>

struct tms tms_buffer;
double ST, UT, st, ut, tt;
int hour, minute;
double second;
long ticks;

#define INIT_TIME \
  ticks = sysconf(_SC_CLK_TCK);\
  times (&tms_buffer);\
  ST = (double) tms_buffer.tms_stime ;\
  UT = (double) tms_buffer.tms_utime

#define FORMAT_TIME_AND_PRINT_TIME \
  times (&tms_buffer);\
  st = ((double) tms_buffer.tms_stime - ST)/ ticks;\
  ut = ((double) tms_buffer.tms_utime - UT)/ ticks;\
  tt = st + ut;\
  hour = (int) (tt / 3600.0);\
  minute = (int) (tt / 60.0) - hour * 60;\
  second =  (tt) - 60 * (int) (tt / 60.0);\
  printf("(%dh %dm %fs)\n", hour, minute, second);

using namespace std;

int Dm, Dn, Dp, Dk, *VecA, *VecC,*VecD,*VecE,**MatD;
double **MatA,**MatB,**MatV;
std::vector <vector <double> > MatY;
std::vector <vector <int> > MatX;

void readValue(char *fn1){
  ifstream fin;
  int i, j, l, v;
  double s;
  char f;
  fin.open(fn1);
  if(!fin){
    cerr << "ERROR: Failed to open file" << endl;
    exit(1);
  }
  fin >> Dm >> Dn >> Dk;
  VecA = new int [Dm];
  VecE = new int [Dm];
  MatX.resize(Dm);
  MatY.resize(Dm);
  for(i = 0; i < Dm; i++){
    fin >> VecA[i];
    for(j = 0; j < VecA[i]; j++){
      fin >> v >> f >> s;
      MatX[i].push_back(v-1);
      MatY[i].push_back(s);
    }
    fin >> VecE[i];
  }
  fin.close();
}

void initValue(){
  int i, j, k;
  double v;
  for(i = 0; i < Dm; i++){
    for(j = 0, v = 0.0; j < VecA[i]; j++){
      v += MatY[i][j]*MatY[i][j];
    }
    for(j = 0, v = 1.0/sqrt(v); j < VecA[i]; j++) MatY[i][j] *= v;
  }//normalize
  MatA = new double *[Dm*(Dm-1)/2];//Similarity matrix
  MatB = new double *[2];
  VecC = new int [Dm*(Dm-1)/2];
  for(i=0; i<Dm*(Dm-1)/2; i++) MatA[i] = new double[3];
  for(i = 0; i < 2; i++){
    MatB[i] = new double [Dn];
    for(j=0;j<Dn;j++) MatB[i][j] = 0;
  }
}

double calSimilarity(int x,int y)
{
  int i,j,k;
  double sum,w;
  for(i=0;i<VecA[x];i++) MatB[0][MatX[x][i]] = MatY[x][i];
  for(i=0;i<VecA[y];i++) MatB[1][MatX[y][i]] = MatY[y][i];
  for(i=0,sum=0.0;i<Dn;i++){
    if(MatB[0][i] != 0 && MatB[1][i] != 0){
      sum += MatB[0][i]*MatB[1][i];
    }
  }
  for(i=0;i<VecA[x];i++) MatB[0][MatX[x][i]] = 0;
  for(i=0;i<VecA[y];i++) MatB[1][MatX[y][i]] = 0;
  return sum;
}

void calValue(){
  int i, j, k=0, l;
  double v,w;
  for(i = 0; i < Dm; i++){
    for(j = i+1; j < Dm; j++){
      MatA[k][0] = i+1; MatA[k][1] = j+1;
      MatA[k][2] = calSimilarity(i,j);//Similarity Vector
      if(k == 0 || MatA[k][2] > v) v = MatA[k][2];
      if(k == 0 || MatA[k][2] < w) w = MatA[k][2];
      k++;
    }
  }
  for(i = 0, v = (Dp-1)/(v-w); i < Dm*(Dm-1)/2; i++) VecC[i] = (int) floor(v*(MatA[i][2]-w));
}

void swapMember(int k, int i, int j)
{
  int h;
  double x;
  h = MatD[k][i]; MatD[k][i] = MatD[k][j]; MatD[k][j] = h; 
  x = MatV[k][i]; MatV[k][i] = MatV[k][j]; MatV[k][j] = x; 
}

void mySort(int k)
{
  int h, i, j;
  double x;
  for(i = 1; i < VecD[k]; i++){
    if(MatV[k][i-1] >= MatV[k][i]) continue;
    swapMember(k, i-1, i); 
    for(j = i-1; j > 0; j--){
      if(MatV[k][j-1] >= MatV[k][j]) break;
      swapMember(k, j-1, j); 
    }
  }
}

void printVector(char *fn1, char *fn2){
  ofstream fout;
  int i, j, k, v;
  VecD = new int [Dp];
  MatD = new int *[Dp];
  MatV = new double *[Dp];
  for(j = 0; j < Dp; j++) VecD[j] = 0; 
  for(i = 0; i < Dm*(Dm-1)/2; i++) VecD[VecC[i]] += 1;
  for(j = 0; j < Dp; j++){ 
    if(VecD[j] == 0) continue;
    MatD[j] = new int [VecD[j]];
    MatV[j] = new double [VecD[j]];
  }
  for(j = 0; j < Dp; j++) VecD[j] = 0; 
  for(i = 0; i < Dm*(Dm-1)/2; i++){ 
    k = VecC[i]; 
    MatD[k][VecD[k]] = i; 
    MatV[k][VecD[k]++] = MatA[i][2]; 
  }
  for(j = 0; j < Dp; j++) if(VecD[j] > 1) mySort(j); 
  fout.open(fn1);
  fout << Dm << " " << Dm*(Dm-1)/2 << endl;
  for(i = Dp-1; i >= 0; i--){ 
    if(VecD[i] == 0) continue;
    for(j = 0; j < VecD[i]; j++){  
      k = MatD[i][j]; 
      fout << MatA[k][0] << " " << MatA[k][1] << " " << MatA[k][2] << endl;
    }
  }
  fout.close();
  fout.open(fn2);
  for(i=0;i<Dm;i++) fout << VecE[i] << endl;
  fout.close();
}


int main(int argc, char **argv){
  int i, j, k;
  double v;
  srand(time(0));
  INIT_TIME;
  Dp = 100000;
  readValue(argv[1]);
  cout << "read OK\n";
  initValue();
  cout << "init OK\n";
  calValue();
  delete [] MatB;
  cout << "cal OK\n";
  printVector(argv[2],argv[3]);
  cout << "print OK\n";
  FORMAT_TIME_AND_PRINT_TIME;
  delete [] VecA;
  delete [] VecD;//Dm
  delete [] VecC;//Dm
  delete [] VecE;
  delete [] MatA;
  delete [] MatV;
  delete [] MatD;
}

ソートの仕方はいろいろありますが、今回はバケットソートを少しいじったものを使いました。

入力
100 784 2
64 128:38 129:254 130:109 156:87 157:252 158:82 184:135 185:241 211:45 212:244 213:150 239:84 240:254 241:63 267:202 268:223 269:11 294:32 295:254 296:216 322:95 323:254 324:195 350:140 351:254 352:77 377:57 378:237 379:205 380:8 405:124 406:255 407:165 433:171 434:254 435:81 460:24 461:232 462:215 488:120 489:254 490:159 516:151 517:254 518:142 544:228 545:254 546:66 571:61 572:251 573:254 574:66 599:141 600:254 601:205 602:3 626:10 627:215 628:254 629:121 654:5 655:198 656:176 657:10 1
・
・
・

出力
2115♯文字数 2235555♯N*(N-1)/2
1070 1392 0.987203
252 810 0.985827
942 1316 0.983981
1147 1983 0.983843
36 279 0.983543
1644 1646 0.982974
1134 1774 0.982413
1802 1821 0.982156
1294 1457 0.982146
・
・
・

続いてこのデータからK-nnアルゴリズムを用いて、各文字に対して類似度の近いk個の文字を得、得られた中で「０」と「１」どちらが多いか調べ、各文字と同じ文字を選択できているか判定します。

結果として、　
「0」成功率　失敗率
「1」成功率　失敗率
という2×2の行列で出力します。

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <vector>
#include <algorithm>
#include<sys/times.h>

struct tms tms_buffer;
double ST, UT, st, ut, tt;
int hour, minute;
double second;
long ticks;

#define INIT_TIME \
  ticks = sysconf(_SC_CLK_TCK);\
  times (&tms_buffer);\
  ST = (double) tms_buffer.tms_stime ;\
  UT = (double) tms_buffer.tms_utime

#define FORMAT_TIME_AND_PRINT_TIME \
  times (&tms_buffer);\
  st = ((double) tms_buffer.tms_stime - ST)/ ticks;\
  ut = ((double) tms_buffer.tms_utime - UT)/ ticks;\
  tt = st + ut;\
  hour = (int) (tt / 3600.0);\
  minute = (int) (tt / 60.0) - hour * 60;\
  second =  (tt) - 60 * (int) (tt / 60.0);\
  printf("(%dh %dm %fs)\n", hour, minute, second);

using namespace std;

int Dm, Dn, Dk, Dl,*VecA, *VecB,*VecD,*VecG,*VecH;
double *VecC,**MatA;
std::vector <vector <int> > MatX;

void readValue(char *fn1, char *fn2){
  ifstream fin;
  int i, j, l, v, f;
  double s;
  fin.open(fn1);
  if(!fin){
    cerr << "ERROR: Failed to open file" << endl;
    exit(1);
  }
  fin >> Dm >> Dn;
  VecA = new int [Dn];
  VecB = new int [Dn];
  VecC = new double [Dn];
  VecG = new int [Dm];
  for(i=0;i<Dl;i++) VecD[i] = 0;
  MatX.resize(Dm);
  for(i = 0; i < Dn; i++){
    fin >> v >> f >> s;
    VecA[i] = v-1;
    VecB[i] = f-1;
    VecC[i] = s; 
  }
  fin.close();
  fin.open(fn2);
  fin >> Dm >> Dl;
  VecD = new int [Dl];
  for(i=0;i<Dm;i++){
    fin >> VecG[i];
    VecD[VecG[i]] += 1;
  }
  fin.close();
}

void initValue(){
  int i, j, k;
  double v;
  VecH = new int [Dl];
  MatA = new double *[Dl];
  for(i=0;i<Dl;i++){
    MatA[i] = new double [Dl];
  }
}

void calValue(){
  int i, j, k, x, y;
  for(i = 0, k = 0; i < Dn; i++){
    x = VecA[i]; y = VecB[i];
    if(MatX[x].size() < Dk){
      MatX[x].push_back(y);
      if(MatX[x].size() == Dk) k+=1;
    }
    if(MatX[y].size() < Dk){
      MatX[y].push_back(x);
      if(MatX[y].size() == Dk) k+=1;
    }
    if(k==Dm) break;
  }
}

void check()
{
  int i,j,k;
  for(i=0;i<Dm;i++){
    for(j=0;j<Dk;j++){
      VecH[VecG[MatX[i][j]]] += 1;
    }
    k = (VecH[0] > VecH[1])? 0:1;
    MatA[VecG[i]][k]+=1;
    for(j=0;j<Dl;j++) VecH[j] = 0;
    if((VecG[i] == 0 && k == 1) || (VecG[i] == 1 && k == 0)){
      printf("%d ",i+1);
      for(j=0;j<Dk;j++) printf("%d ",MatX[i][j]+1);
      printf("\n");
    }
  }
}

void printValue()
{
  int i,j;
  for(i=0;i<Dl;i++){
    for(j=0;j<Dl;j++){
      printf("%f ",MatA[i][j]/VecD[i]);
    }
    printf("\n");
  }
}

int main(int argc, char **argv){
  int i, j, k;
  double v;
  srand(time(0));
  INIT_TIME;
  Dk = atoi(argv[1]);
  readValue(argv[2],argv[3]);
  cout << "read OK\n";
  initValue();
  cout << "init OK\n";
  calValue();
  cout << "cal OK\n";
  check();
  cout << "trans OK\n";
  printValue();
  cout << "print OK\n";
  FORMAT_TIME_AND_PRINT_TIME;
  delete [] VecA;
  delete [] VecB;
  delete [] VecC;//Dm 
  delete [] VecD;
  delete [] VecH;
  delete [] MatA;
}

出力
0.997959 0.002041
0.000881 0.999119
print OK
(0h 0m 2.100000s)

「０」と「１」は全然違うのでほぼ100%の性能です。これが「１」と「７」なんかになるとだいぶ下がってきます。

ちなみに誤判定したものの例を挙げると

確かに誤判定しそうだなぁという感じ。
せっかくなのでパラメトリックな手法と比較してみた方がよかったですね。