データ解析

はじめに

データ解析には、いくつかの標準要素があります。

前処理 — 可能なモデルを同定するために、外れ値、欠損値、平滑化したデータを考える。
まとめ — 全体的な位置、スケール、データの形状を記述する基本的な統計量を計算する。
可視化 — パターンと傾向を確認するために、データをプロットする。
モデリング — 新しい値の予測に適したデータ傾向を詳しく記述する。

データ解析は、以下の 2 つの基本的な目的をもち、上記のことを行います。

正確な予測を行える単純なモデルを用いてデータのパターンを記述する。
モデルを導く変数間の関係を理解する。

この節では、MATLAB^® 環境で基本的なデータ解析を行う方法を説明します。

データの前処理

ライブスクリプトを開く

この例では、解析データを前処理する方法を説明します。

概要

データを適切な MATLAB® コンテナー変数に読み込み、"不適切な" データと "適切な" データを並べ替えることによりデータ解析を開始します。これは、解析のその後の過程で意味のある結果が得られることを保証するための準備の手順です。

データの読み込み

はじめに、count.dat のデータを読み込みます。

load count.dat

24 行 3 列の配列 count は、3 つの交差点それぞれが各列に対応し、ある一日の 1 時間ごとの交通量を行に含みます。

欠損データ

MATLAB NaN (Not a Number: 数値ではない) の値は通常、欠損データを表すために使用されます。NaN の値は欠損データのある変数が構造体 (以下では、3 つすべての交点で一致するインデックスをもつ 24 行 1 列のベクトル) を維持できるようにします。

関数 isnan を用いて、3 つ目の交点におけるデータの NaN 値をチェックします。

c3 = count(:,3); % Data at intersection 3
c3NaNCount = sum(isnan(c3))

c3NaNCount = 0

isnan は、c3 と同じサイズの logical ベクトルを出力します。このベクトルには、データ内の 24 要素のそれぞれについて、NaN の値が存在するか (1)、または存在しない (0) かが記入されています。この場合、logical 値の総和が 0 になるので、データに値 NaN はありません。

NaN の値は「外れ値」の節でデータに導入されます。

外れ値

外れ値は、それ以外のデータとは傾向が著しく異なるデータ値です。外れ値は、測定誤差により発生したものである場合もあれば、データの重要な特徴を表している場合もあります。外れ値を特定し、それらをどのように取り扱うかは、データとその原因を理解してから決めます。

外れ値を特定するための一般的な方法の 1 つは、平均 $μ$ から一定数の標準偏差 $σ$ 分離れた値を探すことです。以下のコードは、 $μ$ と $η = 1, 2$ の場合の $μ + η σ$ のラインと共に、3 つ目の交差点でのデータのヒストグラムをプロットします。

h = histogram(c3,10); % Histogram
N = max(h.Values); % Maximum bin count
mu3 = mean(c3); % Data mean
sigma3 = std(c3); % Data standard deviation

hold on
plot([mu3 mu3],[0 N],'r','LineWidth',2) % Mean
X = repmat(mu3+(1:2)*sigma3,2,1);
Y = repmat([0;N],1,2);
plot(X,Y,'Color',[255 153 51]./255,'LineWidth',2) % Standard deviations
legend('Data','Mean','Stds')
hold off

このプロットから、いくつかのデータが平均を標準偏差の 2 倍を超えて離れていることがわかります。これらのデータを (特徴としてではなく) 誤差として扱う場合は、これらを以下のように値 NaN で置き換えます。

outliers = (c3 - mu3) > 2*sigma3;
c3m = c3; % Copy c3 to c3m
c3m(outliers) = NaN; % Add NaN values

平滑化とフィルター処理

3 つ目の交点でのデータ (「外れ値」で除かれた外れ値をもつ) の時系列プロットは、以下のようになります。

plot(c3m,'o-')
hold on

20 時間での NaN の値は、プロットにおいてギャップとして現れます。この NaN 値の処理は、MATLAB のプロット関数の特徴です。

ノイズのあるデータは、期待値がランダムに変動します。モデルを構築する前に、データの主な特徴を明らかにするためにデータを平滑化することをお勧めします。平滑化する際には、2 つの基本的な仮定をします。

- 予測変数 (時間) と応答 (交通量) の関係が滑らかである。

- 平滑化アルゴリズムによってノイズが減少するため、推定される期待値が改善される。

MATLAB 関数 convn を用いて、単純移動平均平滑化をこのデータに適用します。

span = 3; % Size of the averaging window
window = ones(span,1)/span; 
smoothed_c3m = convn(c3m,window,'same');

h = plot(smoothed_c3m,'ro-');
legend('Data','Smoothed Data')

平滑化の範囲は、変数 span を用いて制御します。平均化の計算では、平滑化ウィンドウがデータに値 NaN を含むときは、値 NaN を返すので、平滑化データのギャップのサイズが増加します。

関数 filter もデータを平滑化するために使用されます。

smoothed2_c3m = filter(window,1,c3m);

delete(h)
plot(smoothed2_c3m,'ro-','DisplayName','Smoothed Data');

平滑化データは、以前のプロットからシフトします。'same' パラメーターを用いた convn は、データと同じ長さでたたみ込みの中央部分を出力します。filter は、データと同じ長さでたたみ込みの初期部分を出力します。それ以外には、アルゴリズムは同じです。

平滑化は、予測の各値において応答値の分布の中心を推定します。このため、多くの近似アルゴリズムの基本的な仮定、"予測の各値での誤差は独立であること" が無効になります。したがって、モデルを "同定" するためには平滑化されたデータを使用しますが、モデルを "近似" するために平滑化されたデータを使用することは避けます。

データのまとめ

ライブスクリプトを開く

この例では、データをまとめる方法を説明します。

概要

多くの MATLAB® 関数によって、全体的な位置、スケール、データサンプルの形状を把握することができます。

MATLAB® の優れた機能の 1 つは、関数によって 1 つのスカラー値だけではなくデータ配列全体を操作できることです。このことは、関数が "ベクトル化" されているといわれます。ベクトル化によって、配列データを用いて効率的に問題を表すことや、ベクトル化された統計関数を用いた効率的な計算が可能になります。

位置の測度

「標準的な」値を見つけることによりデータサンプルの位置を把握します。位置または「中心傾向」の通常の計測は、関数 mean、median および mode により計算されます。

load count.dat
x1 = mean(count)

x1 = 1×3

   32.0000   46.5417   65.5833

x2 = median(count)

x2 = 1×3

   23.5000   36.0000   39.0000

x3 = mode(count)

x3 = 1×3

    11     9     9

他の統計関数と同様に、上記の MATLAB® 関数は変数を列に保持しながら、観測値ごとにデータを行にまとめます。関数は、3 つの交差点それぞれでのデータの位置を 1 回の呼び出しで計算します。

スケールの尺度

データサンプルのスケールまたは「ばらつき」を測定する多くの方法があります。MATLAB® 関数 max、min、std および var はいくつかの一般的な測定値を計算します。

dx1 = max(count)-min(count)

dx1 = 1×3

   107   136   250

dx2 = std(count)

dx2 = 1×3

   25.3703   41.4057   68.0281

dx3 = var(count)

dx3 = 1×3
10³ ×

    0.6437    1.7144    4.6278

他の統計関数と同様に、上記の MATLAB® 関数は変数を列に保持しながら、観測値ごとにデータを行にまとめます。関数は、3 つの交差点それぞれでのデータのスケールを 1 回の呼び出しで計算します。

分布の形状

分布の形状は、分布の位置またはスケールに比べて説明し難いです。MATLAB® 関数 hist のヒストグラムのプロットは、概要を視覚的に表します。

figure
hist(count)
legend('Intersection 1',...
       'Intersection 2',...
       'Intersection 3')

パラメトリックモデルでは、分布の形状の解析的に把握できます。データ平均のパラメーター mu をもつ指数分布は、交通量のデータに対する選択として適切です。

c1 = count(:,1); % Data at intersection 1
[bin_counts,bin_locations] = hist(c1);
bin_width = bin_locations(2) - bin_locations(1);
hist_area = (bin_width)*(sum(bin_counts));

figure
hist(c1)
hold on

mu1 = mean(c1);
exp_pdf = @(t)(1/mu1)*exp(-t/mu1); % Integrates
                                   % to 1
t = 0:150;
y = exp_pdf(t);
plot(t,(hist_area)*y,'r','LineWidth',2)
legend('Distribution','Exponential Fit')

一般のパラメトリックモデルをデータ分布に近似する方法は、この節では取り扱いません。Statistics and Machine Learning Toolbox™ ソフトウェアでは、分布パラメーターの最尤推定値を計算する関数が提供されています。

概要

データのパターンと傾向を可視化するために、多くの種類の MATLAB グラフを利用できます。この節で述べる散布図によって、異なる交差点での交通量データの関係を可視化できます。データ探索ツールによって、グラフ上の個々のデータ点を参照したり、対話形式での操作が行えます。

メモ

この節では、データのまとめのデータ解析を続けます。

2 次元散布図

関数 scatter で作成される 2 次元散布図は、最初の 2 つの交差点における交通量の関係を説明します。

load count.dat
c1 = count(:,1); % Data at intersection 1
c2 = count(:,2); % Data at intersection 2

figure
scatter(c1,c2,'filled')
xlabel('Intersection 1')
ylabel('Intersection 2')

関数 cov で計算される "共分散" は、2 変数間の線形関係の強さを測定します。散布図から最小二乗線に沿ってデータがどれほど密にあるかを評価します。

C12 = cov([c1 c2])

C12 = 2×2
10³ ×

    0.6437    0.9802
    0.9802    1.7144

結果は、対称正方行列として表示されます。(i, j) 番目の位置の要素は、i番目の変数と j番目の変数の共分散です。i番目の対角要素は、i番目の変数の分散です。

共分散は、個々の変数の測定に用いる単位に依存するという不都合があります。共分散の値は、変数の標準偏差で除算することで +1 と -1 の間に正規化できます。関数 corrcoef では "相関係数" を計算します。

R12 = corrcoef([c1 c2])

R12 = 2×2

    1.0000    0.9331
    0.9331    1.0000

r12 = R12(1,2) % Correlation coefficient

r12 = 0.9331

r12sq = r12^2 % Coefficient of determination

r12sq = 0.8707

相関係数の値は正規化されているので、交差点の他の組の値と容易に比較できます。相関係数の二乗である "決定係数" は、最小二乗ラインからの変動を平均からの変動で除算したものです。したがって、決定係数は応答 (この場合、交差点 2 の交通量) における変動の部分であり、散布図から削除されるか、または最小二乗ラインによって統計的に説明されます。

3 次元散布図

関数 scatter3 で作成される 3 次元散布図は、3 つの交差点すべてにおける交通量間の関係を示します。前の手順で作成した変数 c1、c2、c3 を使用します。

figure
c3 = count(:,3); % Data at intersection 3
scatter3(c1,c2,c3,'filled')
xlabel('Intersection 1')
ylabel('Intersection 2')
zlabel('Intersection 3')

関数 eig で共分散行列の固有値を計算することにより、3 次元散布の変数間の線形関係の強さを測定します。

vars = eig(cov([c1 c2 c3]))

vars = 3×1
10³ ×

    0.0442
    0.1118
    6.8300

explained = max(vars)/sum(vars)

explained = 0.9777

固有値は、データの "主成分" の分散です。変数 explained は、データの軸に沿った、第 1 主成分によって説明される変動の割合を測ります。2 次元散布に対する決定係数とは異なり、この尺度では予測変数と応答変数が区別されます。

散布図配列

関数 plotmatrix を使用して、交差点の複数の組間の関係を比較します。

figure
plotmatrix(count)

配列の (i, j) 番目の位置のプロットは、垂直軸上の i 番目と水平軸上の j 番目の変数の散布図です。i 番目の対角位置のプロットは、i番目の変数のヒストグラムです。

グラフ内のデータの探査

ほとんどの MATLAB グラフの観測値は、Figure ツールバーの 2 つのツールを利用してマウスで選択できます。

データカーソル
データのブラシ選択

これらの各ツールはどれも探査モードでの利用となり、そこでグラフ上のデータ点を選択して値を確認したり、特定の観測値を含むワークスペース変数を作成したりできます。データのブラシ選択を使った場合、選択した観測値のコピー、削除または置き換えも可能です。

たとえば、count の第 1 列と第 3 列の散布図を作成します。

load count.dat
scatter(count(:,1),count(:,3))

データカーソルツール

を選択し、右端のデータ点をクリックします。点の x と y の値を表示するデータは、以下にあります。

既定の設定では、データヒントには x 座標、y 座標および z 座標 (3 次元プロットの場合) が表示されます。あるデータ点から他のデータ点にデータをドラッグして新規の値を見るか、データを右クリックしてコンテキストメニューを利用して、データヒントを追加します。MATLAB コードを使用して、データヒントで表示するテキストをカスタマイズすることもできます。

"データのブラシ選択" は、クリックまたはドラッグによって、グラフ上の 1 つまたは複数の観測値を強調表示できるようにした関連機能です。データのブラシ選択モードに入るには、Figure ツールバーで、[データのブラシ選択] ツールの左側をクリックします。ツールアイコンの右側の矢印をクリックすると、観測値のカラーを選択するためのドロップダウンのカラーパレットが開きます。次の Figure は前の Figure と同じ散布図ですが、標準偏差を超えるすべての観測値 ([ツール] 、 [データの統計] GUI により特定) が赤でブラシ選択されています。

scatter(count(:,1),count(:,3))

データの観測値をブラシ選択した後、次の操作を実行できます。

データの観測の削除
データの観測を定数値で置き換え
データの観測を NaN 値で置き換え
データの観測を、コマンドウィンドウにドラッグ、コピー、貼り付け
データの観測をワークスペース変数として保存

たとえば、データのブラシ選択のコンテキストメニューまたは [ツール] 、 [ブラシ選択] 、 [新規変数の作成] オプションを使用して、count13high と呼ばれる新しい変数を作成します。

ワークスペースの新しい変数は、次のようになります。

count13high

count13high =
    61   186
    75   180
   114   257

"リンク付きプロット" または "データリンク" は、データのブラシ選択に関連する機能です。プロットは、プロットが描くワークスペースデータに応答する接続をもつ場合に、リンクしているといわれます。オブジェクトの XData、YData (必要に応じて ZData) に保存された変数のコピーは、これらがリンクしているワークスペース変数が変更または削除されると自動的に更新されます。そのため、これらの変数を表示しているグラフが自動的に更新されます。

変数にプロットをリンクすると、さまざまな表示で特定の観測値をトラックできます。リンクされたプロットのデータ点をブラシ選択する場合、1 つのグラフでブラシ選択すると、同じ変数にリンクしている各グラフ上の同じ観測値が強調表示されます。

データリンクは、変数エディターがワークスペース変数とやりとりするのと同じように、Figure とワークスペース変数の 2 方向の直接のやりとりを確立します。Figure のツールバー上にあるデータリンクツールをアクティブにしてリンクを作成します。このツールをアクティブにすると、次の図に示すリンクプロットのメッセージバー (おそらくタイトルは非表示) がプロットの上部に表示されます。プロットとリンク解除せずに (以下の図に示す) メッセージバーを非表示にできます。この場合、メッセージバーは表示されず、Figure と共に保存されません。

以下の 2 つのグラフは、左のグラフのいくつかの観測値をブラシ選択した後、リンクしたデータの散布図を表します。共通の変数 count は、右の Figure にブラシ選択のマークを置きます。右のグラフは、データのブラシ選択モードではなくても、その変数にリンクしているのでブラシ選択のマークを表示します。

figure
scatter(count(:,1),count(:,2))
xlabel ('count(:,1)')
ylabel ('count(:,2)')
figure
scatter(count(:,3),count(:,2))
xlabel ('count(:,3)')
ylabel ('count(:,2)')

右のプロットは、ブラシ選択された観測値が左のプロットに比べてより直線的な並びであることを示しています。

ブラシ選択されたデータ観測値は、次に示すように、変数エディターにこれらの変数を表示すると、ブラシカラーで強調表示されます。

openvar count

変数エディターにおいて、リンクしたプロットデータの任意の値を変更でき、グラフが編集結果を反映します。変数エディターからデータ観測値をブラシ選択するには、[ブラシ選択ツール] ボタンをクリックします。ブラシ選択した変数がリンクプロットに現在描かれている場合には、ブラシ選択した観測値は変数エディター同様、プロットでも強調表示されます。行列の列である変数をブラシ選択すると、その行内のその他の列もブラシ選択されます。つまり、行ベクトルまたは列ベクトル内の個々の観測値のブラシ選択はできますが、クリックした観測値だけでなく、行列内のすべての列がブラシ選択行で強調表示になります。

データのモデリング

概要
多項式回帰
一般線形回帰

概要

パラメトリックモデルは、データの関係について理解した内容を予測能力をもつ解析ツールへと変換します。トラフィックデータの上昇または下降トレンドに対して、単純に多項式モデルまたは正弦波モデルを選択します。

多項式回帰

関数 polyfit を使用して多項式モデルの係数を推定し、次に、関数 polyval を使用して予測子の任意の値でモデルを評価します。

次のコードは、6 次の多項式モデルで 3 つ目の交点でのトラフィックデータを近似します。

load count.dat
c3 = count(:,3); % Data at intersection 3 
tdata = (1:24)'; 
p_coeffs = polyfit(tdata,c3,6);

figure 
plot(c3,'o-') 
hold on 
tfit = (1:0.01:24)'; 
yfit = polyval(p_coeffs,tfit); 
plot(tfit,yfit,'r-','LineWidth',2)
legend('Data','Polynomial Fit','Location','NW')

このモデルは、単純ながら上下のトレンドに追従する利点があります。ただし、特にデータの両端など、予測能力の精度には疑問があります。

一般線形回帰

データには、12 時間の周期があり、7 時付近にピークがあることを仮定すると、次の形の正弦波モデルで近似するのが適切です。

$y = a + b \cos ((2 π / 12) (t - 7))$

係数 a と b は直線的に表示されます。MATLAB® mldivide (バックスラッシュ) 演算子を使用して一般的な線形モデルへ近似します。

load count.dat
c3 = count(:,3); % Data at intersection 3 
tdata = (1:24)'; 
X = [ones(size(tdata)) cos((2*pi/12)*(tdata-7))];
s_coeffs = X\c3;

figure
plot(c3,'o-')
hold on
tfit = (1:0.01:24)';
yfit = [ones(size(tfit)) cos((2*pi/12)*(tfit-7))]*s_coeffs; 
plot(tfit,yfit,'r-','LineWidth',2)
legend('Data','Sinusoidal Fit','Location','NW')

関数 lscov を使用して、係数の推定標準誤差や平均二乗誤差など、近似の統計値を計算します。

[s_coeffs,stdx,mse] = lscov(X,c3)

s_coeffs = 2×1

   65.5833
   73.2819

stdx = 2×1

    8.9185
   12.6127

mse = 1.9090e+03

データの 12 時間間隔の仮定を、関数 fft を使用して計算した "ピリオドグラム" によって確認します。

Fs = 1; % Sample frequency (per hour)
n = length(c3); % Window length
Y = fft(c3); % DFT of data
f = (0:n-1)*(Fs/n); % Frequency range
P = Y.*conj(Y)/n; % Power of the DFT

figure
plot(f,P)
xlabel('Frequency')
ylabel('Power')