トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索   ヘルプ   最終更新のRSS

MATLAB Note/統計

Last-modified: 2013-03-07 (木) 05:28:34
Top / MATLAB Note / 統計

教師あり学習 †

↑

決定木 †

  • まずは入力値データ（独立変数）と、結果データ（従属変数）を宣言します。
    %独立変数（天気、気温、湿度、風の強さ）
meas = [
   1 29 85 1
   1 27 90 2
   2 28 78 1
   3 21 96 1
   3 20 80 1
   3 18 70 2
   2 18 65 2
   1 22 95 1
   1 21 70 1
   3 24 80 1
   1 24 70 2
   2 22 90 2
   2 27 75 1
   3 22 80 2
]
%従属変数（ゴルフをするか、しないか）
species = [
   '×'
   '×'
   '○'
   '○'
   '○'
   '×'
   '○'
   '×'
   '○'
   '○'
   '○'
   '○'
   '○'
   '×'    
]
  • 決定木の表示は、以下のようにします（上のコードに続けて実行します）。
    t = classregtree(meas, species,'names',{'天気' '気温' '湿度' '風'})
view(t);
  • 結果は以下のようになります。
    kettei_1.JPG
  • 上記の結果は入力値の数が足りないため、十分な結果とはいえません。そこで、measとspeciesの値を増やして、もう一度決定木をプロットします（上のコードに続けて実行します）。
    meas2 = [meas ; meas];                     %measの値を2度繰り返し与える
species2 = [species ; species];            %speciesの値を2度繰り返し与える
t2 = classregtree(meas2, species2,'names',{'天気' '気温' '湿度' '風'})
view(t2);
    kettei_2.JPG 
    meas3 = [meas ; meas ; meas];              %measの値を3度繰り返し与える
species3 = [species ; species ; species];  %speciesの値を3度繰り返し与える
t3 = classregtree(meas3, species3,'names',{'天気' '気温' '湿度' '風'})
view(t3);
    kettei_3.JPG 
    meas5 = [meas ; meas ; meas ; meas ; meas];
species5 = [species ; species ; species ; species ; species];
t5 = classregtree(meas5, species5,'names',{'天気' '気温' '湿度' '風'})
view(t5);
    kettei_4.JPG
  • 入力値を3回以上繰り返し与えると、決定木が収束しました。
  • 「湿度が高い晴れの日」と「風が強い雨の日」はゴルフをしない（×）、という結果が求められました。
  • 決定木の評価は 関数 TEST EVAL ISBRANCH CUTPOINT CUTVAR CUTCATEGORIES CUTTYPE NODEERR NODEPROB NODESIZE PARENT PRUNE RISK などを参照してください（コマンドラインで呼び出すには classregtree/関数名 ）。
↑

サポートベクターマシン - Support Vector Machine Toolbox*1 †

  • 【注意】
    • いまのところ、32ビット Windows環境でしか動きません。qp.mexw32 を以下の作業フォルダ内に入れておいてください。
  • 導入
    • リンク先のZIPファイルをダウンロードして、解凍してできた「svm」フォルダを、MATLAB作業フォルダの適当なところにおきます。
    • フォルダのパスを指定します。「svm」フォルダの一階層上のフォルダで、以下のコードを実行してください。 
      %フォルダの実行パスをMATLABサーチパスに加える
addpath 'svm';
      • 「svm」フォルダを右クリック → パスに追加 → 選択フォルダ としてもＯＫです。
      • README には MATLAB のサーチパスに加えるように、とあります。これでも構いません。
  • 実行
    • 以下は全て、「svm」フォルダの一階層上のフォルダで実行しています。
    • 学習データと教師データを指定
      • ここでは、とりあえず、SVM TOOLBOX v0.55 のサンプルデータ iris_inputdata.txt を使ってみます。 
        % 学習データを指定
data = load('iris_inputdata.txt');    %データセットを読み込む
x = data(:, 1:2)                      %学習データ（入力ベクトル）
y = data(:, 3)                        %正解データ（ターゲットベクトル）
% データを表示
disp([x y])
      • 左から２列が入力ベクトル、３列目がターゲットベクトル（1 か -1 か）です。
    • カーネルを指定
      • カーネルはSVMの分類学習のルールのようなものです。このToolboxが提供しているカーネルの詳細は、 
        help svkernel
        で確認できます。なお、存在しない名前を指定した場合、「liner（線形カーネル）」がセットされます。
      • ここでは、スプラインカーネルを指定しておきます。 
        kernel = 'spline';
    • 境界を推定
      • 続いて、境界の推定学習を行ないます。 
        [nsv, alpha, b0] = svc(x, y, kernel)       % x : 入力ベクトル y : ターゲットベクトル
                                           %　詳細は >> help svc を参照
[h] = svcplot(x, y, kernel, alpha, b0)     %結果をプロット
      • 以下のようにプロットされます。
        svm_spline.png
      • 赤丸と青丸がターゲットベクトルの種類、縦軸と横軸が入力ベクトルの１列目と２列目に対応します。
      • SVM TOOLBOX v0.55 の結果とは少し境界が異なりますが、うまくいっているようです。
    • とりあえず、以下のカーネルが動きました（他は謎のエラーが出ます）。
      • erfb カーネルと anova カーネル
        svm_erfb.png
      • liner カーネル
        svm_liner.png
  • どちらも、おそらく飛び地になっている赤が原因で、うまく分類できていません。spline カーネルを使うのがよさそうです。
  • 評価
    • 新しい入力を SVM の境界と照らし合わせて、どちらのターゲットに分類されるかを確認するには、 svcoutput を使います。
    • ここでは、簡単な学習データと教師データを与えて境界を推定し、新たに与えたテストデータがどちらに分類されるかを調べます。 
      %フォルダの実行パスをMATLABサーチパスに加える
addpath 'svm';

% 学習データを指定
x = [1 1 ; 2 2 ; 3 3 ; 4 4];
% 教師データを指定
y = [1   ; 1   ;-1   ;-1  ];
disp([x y]);

% カーネルを指定
kernel = 'spline'; 

% 境界を推定
[nsv, alpha, b0] = svc(x, y, kernel);
% 結果をプロット
[h] = svcplot(x, y, kernel, alpha, b0);

% テストデータを与えて評価
tstX = [1 0 ; 4 3];
tstOUTPUT = svcoutput(x, y, tstX, kernel, alpha, b0);
disp([tstX tstOUTPUT]);
  • 結果は以下のようになります。
    svm_output.png
  • 左上の数字が学習データと教師データ、左下の数字がテストデータとその分類結果です。
  • [1 0] がターゲット 1 （青丸）、[4 3] がターゲット -1 （赤丸）になっていますから、正しく分類できているといえます。
  • 注意事項
    • README によれば、上記のプログラムは MATLAB で書かれていることから、処理速度は遅めのようです。
      • MATLAB のソースは勉強用で、実際のクラスタリングを行なうときは、パッケージに同梱されている C プログラム版をコンパイルして使用することが推奨されています。
      • C 版については、こちらでは動かしていないので、よく分かっていません。申し訳ないです。
    • 入力ベクトルが３次元以上になっても、（svcplot は使えないですが）おそらく上と同じ手続でうまくいくはずです。
    • これがソフトマージン法（多少のエラーを許容して境界を引く）なのか、ハードマージン法（許容しない）なのかは、確認できていません。もしもハードマージンであれば、データの種類によってはエラーが出るかもしれませんが、パッケージ内にソフトマージンの関数があるので、対処可能だと思います。
↑

サポートベクターマシン - MATLAB SVM TOOLBOX v0.55 †

  • SVM TOOLBOX v0.55
    • 【注】これはデモプログラムがあり、機能が分かりやすいかと思ったのですが、MATLAB Support Vector Machine Toolbox の方が完成度は高く、より分かりやすかったため、こちらは中断しています。
  • 導入
    • 解凍してできた「svm_v0.55」フォルダを、MATLAB作業フォルダの適当なところにおきます。
    • このパッケージはC++プログラムも使用しているので、C++のコンパイラがない場合はインストールしておく必要があります。
      • Windowsなら、Borland C++Compilerなどをダウンロードしてインストールしてください。
    • インストールしたC++のコンパイラをMATLABで使えるようにします。以下のコードを実行してください。 
      mex -setup
      • メッセージが表示されますので、[Y]>[(C++Compilerの番号)]>[Y]と入力してください。
  • 実行
    • デモプログラムを実行してみます。「svm_v0.55」フォルダ直下で、以下のコードを実行してください。 
      addpath 'svm_v0.55';   %フォルダの実行パスをMATLABサーチパスに加える
compilemex;            %@smosvctutorフォルダ内のC++（cpp）ファイルをコンパイル
demo                   %demo.mを実行
    • 結果は以下のようになります。デモの結果をそのまま転載したものです。
      svm_demo.jpg
    • demo.m の中身を編集して、SVM TOOLBOXのアルゴリズムを勉強できるようになっています。
↑

教師なし学習 †

↑

主成分分析 †

  • 同等の機能を持った複数の関数が提供されているようです。
    • pcacov (Statistics Toolbox)
    • princomp (Statistics Toolbox)
    • prepca (Neural Network Toolbox) － 実行例
      • prepcaは古い関数です。最新版では、同様の機能を持ったprocesspcaの使用が推奨されています。
      • エラーが出たとき、ここが参考になるかも
  • processpca による方法
    • ここでは、Neural Network Toolboxが提供している主成分分析の手法を紹介します（参考）。
%データをロード
load choles_all
%  データ choles_all の内訳
%   p : 264人の患者の血液サンプルのスペクトルの測定結果（21個のスペクトルの波長）
%   t : 血漿分離に基づくhdl, ldl, vldlコレステロールレベルを測定したデータ → こちらは使いません。
size(p)           %変換前のデータサイズをチェック（21次元264個のデータ）

%主成分分析を実行
pn = prestd(p);                   %平均が0で標準偏差が1になるようにデータを前処理 
ptrans = processpca(pn, 0.001);   %データセット内の99.9%の変動にあたる主成分を取得
size(ptrans)      %変換されたデータサイズをチェック（4次元264個のデータ）
  • ベクトル p が 21 次元から 4 次元に縮んでいます。主成分分析によってデータサイズを減らすことができました。
  • ここで使用したprestdも古い関数のようで、最新版ではmapstdの使用が推奨されています。
    • prestdをそのままmapstdに置き換えるとエラーが出るため、ここではとりあえずそのままにしています。
  • princomp による方法
    • Statistics toolbox による主成分分析（参考） 
      function [ tokuten ] = PCA( inputdata, newinputdata )

%主成分分析の前に，平均が0，分散が1になるように正規化
N = size(inputdata, 1);
m = mean(inputdata);
s = std(inputdata);
inputdata2 = (inputdata - repmat(m, N, 1)) ./ repmat(s, N, 1)

%主成分分析
[w, Z, s2, t2] = princomp(inputdata2)

%第2主成分までを表示する
Z(:, 1 : 2)

%寄与率を第2主成分まで表示する
ruiseki = s2 / sum(s2);
ruiseki(1 : 2)

tokuten = [];
%追加データの主成分得点を出力するとき
if nargin == 2
    tokuten = newinputdata * w;
end
↑

階層型クラスタリング †

  • Statistics Toolbox が提供している機能です。
  • 参考(関数)　参考(詳細)
  • パラメータの詳細は 関数 PDIST LINKAGE DENDROGRAM MANOVACLUSTER のヘルプを参照してください。
  • 例えば、次のような動物の特徴データを階層構造で分類したいとします*2。
    animal.JPG
  • まずは入力値と、入力値のラベル（それぞれのデータ名）を宣言します。
     X = [
     1 0 0 1 0 0 0 0 1 0 0 1 0
     1 0 0 1 0 0 0 0 1 0 0 0 0
     1 0 0 1 0 0 0 0 1 0 0 1 1
     1 0 0 1 0 0 0 0 1 0 0 1 1
     1 0 0 1 0 0 0 0 1 1 0 1 0
     1 0 0 1 0 0 0 0 1 1 0 1 0
     0 1 0 1 0 0 0 0 1 1 0 0 0
     0 1 0 0 1 1 0 0 0 0 1 0 0
     0 1 0 0 1 1 0 0 0 0 1 0 0
     0 1 0 0 1 1 0 1 0 1 1 0 0
     1 0 0 0 1 1 0 0 0 1 0 0 0
     0 0 1 0 1 1 0 0 0 1 1 0 0
     0 0 1 0 1 1 0 1 0 1 1 0 0
     0 0 1 0 1 1 1 1 0 0 1 0 0
     0 0 1 0 1 1 1 1 0 0 1 0 0
     0 0 1 0 1 1 1 0 0 0 0 0 0
 ]
 Xlabel = strvcat('ハト','メンドリ','アヒル','ガチョウ','フクロウ',...
    'タカ','ワシ','キツネ','イヌ','オオカミ','ネコ','トラ','ライオン',...
    'ウマ','シマウマ','ウシ')                   %STRVCAT 文字列を垂直に結合
  • 階層型クラスタリングは、以下のようにします（上のコードに続けて実行します）。
     Y= pdist(X,'cityblock');    %PDIST 観測間の距離
 Z= linkage(Y,'average');    %LINKAGE 重みのない平均距離による階層構造クラスタツリーの作成
 [H,T] = dendrogram(Z,'colorthreshold','default','LABELS',Xlabel);
  • 結果は以下のようになります。
    dendro.jpg
  • クラスタリングの結果は、変数Zに格納されています。 
    Z                           %変数Zの中身を表示

Z =
   14.0000   15.0000         0
    8.0000    9.0000         0
    5.0000    6.0000         0
    3.0000    4.0000         0
   12.0000   13.0000    1.0000
    1.0000    2.0000    1.0000
   20.0000   22.0000    1.5000
   19.0000   23.0000    1.7500
   16.0000   17.0000    2.0000
   10.0000   18.0000    2.0000
   21.0000   25.0000    2.8333
   26.0000   27.0000    3.8000
    7.0000   24.0000    3.8333
   11.0000   28.0000    4.3750
   29.0000   30.0000    8.7619
  • データの意味は以下の通りです。
    • Z の 1 行目 = 「データ番号14（ウマ） と 15（シマウマ）の距離が 0 」
      • データ番号 14＋15 の親ノードの番号を、16 + 1 の 17 番とします。
    • Z の 2 行目 = 「データ番号8（キツネ） と 9（イヌ）の距離が 0 」
      • データ番号 8＋9 の親ノードの番号を、16 + 2 の 18 番とします。
    • Z の 9 行目 = 「データ番号16（ウシ） と 17（ウマ＋シマウマ）の距離が 2 」
  • Z の 1 列目 = リーフの左側、 2 列目 = リーフの右側ですから、それぞれに 0 1 をつけてやれば、「ウマ = 000000」「シマウマ = 000001」などといった表記ができそうです。
    • MATLAB の関数で、すでにこれを実現しているものがあるかもしれません（まだ調べていません）。
  • 「>>help linkage」で詳細な解説を読めます。
↑

k-means法 †

  • Statistics Toolbox が提供している機能です。
  • K-平均法、c-平均法とも呼ばれます。 → ウィキペディア「K平均法」
  • 参考
  • 教師なし、非階層的なクラスタリングアルゴリズムの一種です。入力データの集合をk個のクラスタに分割します。
  • 「>>help kmeans」を実行すると、日本語のヘルプを確認できます。
  • 例えば、以下のような分布をとる入力値を４クラスタに分類したいとします。
    kmeans_1.jpg
  • 上図の入力値は以下のようにして作成しました（平均の異なる４つの正規分布）。
     %正規分布の平均値（分布の中心、最も出やすい値）を決める
 X1nd = [zeros(500,1)+2 zeros(500,1)+2];   %正規分布1の平均値 （2,2）、500行2列
 X2nd = [zeros(500,1)+4 zeros(500,1)+0];   %正規分布2の平均値 （4,0）
 X3nd = [zeros(500,1)-2 zeros(500,1)+2];   %正規分布3の平均値（-2,2）
 X4nd = [zeros(500,1)+2 zeros(500,1)+4];   %正規分布4の平均値 （2,4）
 %データ行列Xを生成（1行1入力、2000行2列のデータ）
 X = [randn(500,2) + X1nd ;...             %500行2列の正規分布を生成して平均値を修正
     randn(500,2) + X2nd ;...
     randn(500,2) + X3nd ;...
     randn(500,2) + X4nd ];
 %プロット
 plot(X(:,1),X(:,2),'.k');
  • k-means法によるクラスタリングは、以下のようにします（上のコードに続けて実行します）。
     clusterNum = 4;                           %クラスタの総数（４クラスに分類）
 %クラスタリングを実行する
 [cidx, ctrs] = kmeans(X, clusterNum);     %k-means法で入力値の集合XをclusterNum個のクラスに分類
 %cidx : 分類後の各入力値のクラスタ番号（2000行1列）
 %ctrs : クラスタの中心の座標（4行2列）
 %結果をプロット
 hold on;
 plot(X(cidx==1,1),X(cidx==1,2),'.',...    %(cidx=1の各行の1列目),(2列目),点でプロット
     X(cidx==2,1),X(cidx==2,2),'.',...
     X(cidx==3,1),X(cidx==3,2),'.',...
     X(cidx==4,1),X(cidx==4,2),'.');
 %上のグラフにクラスタ中心点を追記（ダイヤ、黒枠、緑で塗りつぶし、大きさ8）
 plot(ctrs(:,1),ctrs(:,2),'d','MarkerEdgeColor','k',...
                        'MarkerFaceColor','g',...
                        'MarkerSize',8);
 hold off;
  • 結果は以下のようになります。各クラスタの中心点（緑のダイヤ）が、正規分布の平均値を表現できていることが分かります。
    kmeans_2.jpg
  • クラスタの数を５にするには、以下のようにします。 
     clusterNum = 5;
 [cidx, ctrs] = kmeans(X, clusterNum);
 hold on;
 plot(X(cidx==1,1),X(cidx==1,2),'.',...
     X(cidx==2,1),X(cidx==2,2),'.',...
     X(cidx==3,1),X(cidx==3,2),'.',...
     X(cidx==4,1),X(cidx==4,2),'.',...
     X(cidx==5,1),X(cidx==5,2),'.');
 plot(ctrs(:,1),ctrs(:,2),'d','MarkerEdgeColor','k',...
                        'MarkerFaceColor','g',...
                        'MarkerSize',8);
 hold off;
  • 結果は以下のようになります。
    kmeans_3.jpg
  • 関数kmeansの引数を指定して、クラスタリングのパラメータを設定できます（ >>help kmeans を参照）。
  • 【注意】 kmeans法の結果は、実行するたびに変わります。
  • 適正なクラスタ数を情報量基準に基づいて自動的に推定する x-means法 も提案されています。
↑

自己組織化マップ †

↑

ガウス混合分布モデル †

↑

線形モデル †

↑

確率分布 †

↑

正規分布の検定 †

↑

回帰分析 †

  • 目的変数（従属変数）と説明変数（独立変数）の間の関係を表す式を求める。
    • 説明変数が1つ（回帰分析）
      • 目的変数をy、説明変数をxとしたとき、誤差が最小になるようにデータセットを説明可能な関数 y = ax + b を求めればよい。
      • a（傾き）を、回帰係数 と呼びます。
    • 説明変数が複数（重回帰分析）
  • サンプルデータとして、統計解析ソフト R のサンプルデータ、cars を使います。
    • cars_R.csv
    • データプロット 
      data = csvread( 'cars_R.csv' );
plot( data(:,1), data(:,2), '*' );
title('cars'); xlabel('車の速度'); ylabel('ブレーキをかけた後停止するまでの距離');
      reg_001.png
  • 「ブレーキをかけた後停止するまでの距離」を目的変数、「車の速度」を説明変数として、回帰分析をしてみます。
y = data(:,2);
one = ones( length(data), 1 );
x = [one data(:,1)];
coefficient = regress(y, x)
  • 説明変数のベクトルには、単位行列を追加して 2×(データ数) のベクトルにする必要があります*3。この処理をしないと、原点を通過する回帰直線の傾きのみが出力されます。
  • 結果 
    coefficient =

  -17.5791
    3.9324
  • この結果から求まった回帰式（回帰モデル）は以下の通り。
a math image
  • 回帰直線を引いてみます。 
    hold on;
xdata = 0 : 0.01 : 25;
ydata = coefficient(1) + coefficient(2) * xdata;
plot(xdata, ydata);
hold off;
    reg_002.png
↑

一般化線形モデル †

↑

その他の統計モデル †

↑

隠れマルコフモデル †

  • Statistics Toolbox が提供している機能です。 ※公式で英語の解説、helpコマンドで日本語の解説を読めます。 → 実行例
  • hmmdecode －系列の後方の状態確率を計算
  • hmmestimate －HMM に与えられた状態情報に対するパラメータ推定
  • hmmgenerate －隠れマルコフモデルに対する系列の生成
  • hmmtrain －HMM に対するモデルパラメータの最尤推定量
  • hmmviterbi －系列に対して最も起こりうる状態パスを計算
↑

ベイジアンネットワーク †

*1 サポートベクターマシンについては Support Vector Machine って，なに？（産総研 青葉先生）の解説が面白い！
*2 Teuvo Kohonen (1995).『Self-Organizing Maps』, Springer. (徳高平蔵・堀尾恵一・堀尾正昭・大薮又茂・藤村喜久郎訳(2003). 自己組織化マップ.『シュプリンガー・フェアラーク東京』.).
*3 参考:MATLAB 入門３ MATLAB 基本演算（横浜国際大 永井先生）