概要
MATLAB のテーブルデータを整理していると「X と Y の組み合わせごとに Z の平均とりたいなぁ」と思うこと，ありますよね？これを
meanZ = reduceTableByAxes(T, ["X", "Y"], @mean, "Z")
の一行で実現する関数を実装します．
はじめに
東京農工大学 MATLAB Student Ambassador の松本です．
MATLAB の table は，条件や設定値と，それに対応するデータをまとめて保存するのに便利なデータ型です．一方で，保存したデータを解析するときには，
- 特定の条件を満たす行だけを取り出す
- 複数の条件の組み合わせごとにデータを分割する
- 各条件について平均値などを計算する
といった処理が頻繁に必要になります．これらの処理は論理インデックスや for 文を使って実装できますが，条件が増えたり変更されたりするたびにコードを書き直すのが面倒な場合もあります...
本記事では，table を条件ごとに整理するために実装した3つの関数を実装します！
-
filterTable：指定した条件を満たす行を取り出す -
splitTableByAxes：指定した変数を軸に，その組み合わせごとに table を分割する -
reduceTableByAxes：指定した変数を軸に，その組み合わせごとに集約する（＝平均値など計算する）
コードは以下のリポジトリで公開しています．
また，以下のURLから MATLAB Online 上でコードを実行できます．
以下ではまず，reduceTableByAxes の使用例を示した後，各関数の仕様と実装を順に見ていきます．
例：多項式フィッティングと過適合
例として，多項式フィッティングの結果をまとめたテーブル T を考えます（関数 dummyData は記事末尾に記載しています）．
T = dummyData();
このテーブル T には，図中の train = true のデータに対し，次数 N = 1, 2, ..., 10 の多項式をフィッティングした結果が格納されています．ここで，y_observe はマーカーに示す観測値，y_predict は曲線に示す多項式による予測値，error は両者の絶対誤差です．
N = 1（直線近似）は不十分である一方，N = 10はtrain = trueのデータに過適合して，曲線が暴れています．この中ではN = 3が最も良好な結果のようです．
過適合の様子を確認するために，train = true と train = false のそれぞれについて，次数 N に対する誤差 error の平均値をプロットしましょう．
本記事で実装する reduceTableByAxes 関数を使えば，3 行のコード（＋見栄えを整えるための数行）でこの可視化を実現できます．
[averageError, axisValues] = reduceTableByAxes(T, ["N", "train"], @mean, "error");
bar(axisValues.N, averageError)
legend("train=" + string(axisValues.train), Location="northwest")
% 以下は見栄え調整
yscale("log")
xlabel("Polynomial degree N")
ylabel("Mean absolute error")
grid on
train = true のデータに対する誤差は，次数を増やすにつれて概ね減少する一方，train = false のデータに対する誤差は，次数が大きくなりすぎると増加しています．
実装する関数の仕様
filterTable
Tfilt = filterTable(T, filterVar1, filterValue1, filterVar2, filterValue2, ...)
処理対象の table T から，指定した条件をすべて満たす行を抽出します．条件として使用する変数名 filterVar と，その値 filterValue の組は必要な数だけ続けて指定できます．
例えば，T.N == 3 かつ T.train == false の行は，次のように取り出せます．
Tfilt = filterTable(T, "N", 3, "train", false)
Tfilt =
9×6 table
x y_observe y_predict error train N
____ _________ _________ ________ _____ _
-1.8 -1.8777 -2.4668 0.58917 false 3
-1.6 -1.1235 -1.1061 0.017482 false 3
-0.8 0.9204 0.93454 0.014148 false 3
-0.2 0.56793 0.27946 0.28847 false 3
0.8 -1.1608 -1.0122 0.14865 false 3
1.2 -1.2616 -0.44289 0.81873 false 3
1.4 -0.43702 0.27261 0.70963 false 3
1.6 1.2486 1.3414 0.092848 false 3
2 3.9175 4.737 0.81955 false 3
splitTableByAxes
[Tsplit, axisValues] = splitTableByAxes(T, axisVars)
table T を，axisVars に指定した変数の値の組み合わせごとに分割します．axisVars には，配列の軸として使用する変数名を string 配列として指定します．出力 Tsplit には，分割された table を格納した cell 配列が返されます．
例えば，以下のコードはテーブル T を N および train の二つの軸について分割します．
[Tsplit, axisValues] = splitTableByAxes(T, ["N", "train"])
Tsplit =
10×2 の cell 配列
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
{12×6 table} {9×6 table}
axisValues =
フィールドをもつ struct:
N: [10×1 double]
train: [2×1 logical]
このとき，Tsplit の 1次元目は N，2次元目が train に対応しており，その中には T のうち対応する条件を満たす行を抽出したものが格納されています．axisValues には，各次元がどの値に対応するかが以下の形式で格納されます．
axisValues.N = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
axisValues.train = [false, true]
Tsplit{3,1} は filterTable(T, "N", 3, "train", false) に一致します．
reduceTableByAxes
[reducedValues, axisValues] = reduceTableByAxes(T, axisVars, reduceFcn, dataVar)
table T を axisVars に指定した変数の値の組み合わせごとに分割し，各グループの変数 dataVar に関数 reduceFcn を適用します．出力 reducedValues には各条件の組み合わせに対する reduceFcn の出力を並べた配列になります．splitTableByAxes と同様に，その各次元は axisVars の各要素に対応します．
例えば，N と train の二つの軸について error の平均値を求めるには，次のように記述します．
[averageError, axisValues] = ...
reduceTableByAxes(T, ["N", "train"], @mean, "error")
averageError =
1.2601 1.1891
1.2829 1.2110
0.1521 0.3887
0.1265 0.3734
0.1197 0.4119
0.1193 0.4957
0.1019 0.9103
0.1026 2.5052
0.1070 16.1220
0.0816 169.5976
各関数の実装
以下に各関数の完成版を示します．詳しい解説は折りたたみ内に記載しています．
filterTable の実装
完成版
function [Tfilt, idx] = filterTable(T, filterVar, filterValue)
arguments
T table
end
arguments (Repeating) % 繰り返し引数
filterVar (1,1) string % 変数名なので 1x1 の文字列であるべき
filterValue
end
% 条件を満たす行の検索
idx = true(height(T), 1);
for k = 1:numel(filterVar)
idx = idx & ismember(T.(filterVar{k}), filterValue{k});
end
% 条件を満たす行の抽出
Tfilt = T(idx, :);
end
コードの解説はこちら
arguments (Repeating) で条件を繰り返し受け取る
filterTable では，"N", 3 や "train", false のような「変数名と値の組」を，必要な数だけ受け取る必要があります．このような入力には，arguments (Repeating) が使えます．
function [Tfilt, idx] = filterTable(T, filterVar, filterValue)
arguments
T table
end
arguments (Repeating) % 繰り返し引数
filterVar (1,1) string % 変数名なので 1x1 の文字列であるべき
filterValue
end
繰り返し部分の引数は cell 配列として扱われます．例えば，
filterTable(T, "N", 3, "train", false)
の場合，第2引数以降の繰り返し部分は以下のように格納されます．
filterVar = {"N", "train"};
filterValue = {3, false};
条件を論理インデックスへ追加
続いて，各条件を満たす行を示す論理インデックスを用意し，行を抽出します．まず全要素が true の配列を用意し，次いで各条件を満たす行との論理積（AND）を取っています．
% 条件を満たす行の検索
idx = true(height(T), 1);
for k = 1:numel(filterVar)
idx = idx & ismember(T.(filterVar{k}), filterValue{k});
end
% 条件を満たす行の抽出
Tfilt = T(idx, :);
end
T.(filterVar{k}) は，文字列 filterVar{k} に対応する変数を取り出す操作です．例えば filterVar{k} が "N" なら，T.N と同じ意味になります．
比較には == ではなく ismember を使うことで，複数の値の OR 検索を可能にします．
例えば以下の例では，N が2，3，4のいずれかであり，かつ train が false の行が抽出されます．
Tfilt = filterTable(T, "N", [2 3 4], "train", false);
splitTableByAxes の実装
完成版
function [Tsplit, axisValues] = splitTableByAxes(T, axisVars)
arguments
T table
axisVars (1,:) string
end
% 各軸に含まれる値を取得
numAxes = numel(axisVars);
axisSizes = zeros(1, numAxes);
axisValueList = cell(1, numAxes);
axisValues = struct;
for axisIdx = 1:numAxes
axisVar = axisVars(axisIdx);
values = unique(T.(axisVar), "stable");
axisValues.(axisVar) = values;
axisValueList{axisIdx} = values;
axisSizes(axisIdx) = numel(values);
end
% 条件の全組み合わせを作成
axisCombinations = combinations(axisValueList{:});
axisCombinations.Properties.VariableNames = cellstr(axisVars);
% 各条件を満たすテーブルを抽出
Tsplit = cell(height(axisCombinations), 1);
for combinationIdx = 1:height(axisCombinations)
filterArgs = [
num2cell(axisVars)
table2cell(axisCombinations(combinationIdx, :))
];
Tsplit{combinationIdx} = ...
filterTable(T, filterArgs{:});
end
% セル配列の各次元を軸に対応させる
Tsplit = reshape(Tsplit, [fliplr(axisSizes), 1]); % 多次元配列に変形
Tsplit = permute(Tsplit, [numAxes:-1:1, numAxes + 1]); % 軸の順番を逆転
end
コードの解説はこちら
多重の for 文を使う実装
まずは軸を N と train の二つに限定した，素朴な実装から考えてみます．N と train に含まれる値を unique を使って洗い出し，先ほど作った filterTable を呼び出せばよいので，次のように書けます．
nValues = unique(T.N, "stable"); % [1,2,3,4,5,6,7,8,9,10]
trainValues = unique(T.train, "stable"); % [false, true]
Tsplit = cell(numel(nValues), numel(trainValues));
for nIdx = 1:numel(nValues)
for trainIdx = 1:numel(trainValues)
Tsplit{nIdx, trainIdx} = filterTable(T, "N", nValues(nIdx), "train", trainValues(trainIdx));
end
end
ただし，軸が増えるたびに for ループも増やさなければなりません．splitTableByAxes では，この「全組み合わせを列挙して filterTable を呼ぶ」という処理を，軸の数によらず同じコードで実行できるようにしています．
各軸に含まれる値を集める
まず，axisVars で指定された各変数について，テーブル内に存在する値を unique で取得し，これを axisValueList に格納します．
axisSizes には，各軸に含まれる値の個数を保存します．これは，最後に cell 配列を多次元配列へ変形するときに使用します．splitTableByAxes(T, ["N", "train"]); の場合，axisSizes = [10, 2] となります．出力用の axisValues もこの段階で作成します．
% 各軸に含まれる値を取得
numAxes = numel(axisVars);
axisSizes = zeros(1, numAxes);
axisValueList = cell(1, numAxes);
axisValues = struct;
for axisIdx = 1:numAxes
axisVar = axisVars(axisIdx);
values = unique(T.(axisVar), "stable");
axisValues.(axisVar) = values;
axisValueList{axisIdx} = values;
axisSizes(axisIdx) = numel(values);
end
combinations で条件の全組み合わせを作る
各軸の値が得られたので，次は全条件の組み合わせを列挙します．
% 条件の全組み合わせを作成
axisCombinations = combinations(axisValueList{:});
axisCombinations.Properties.VariableNames = cellstr(axisVars);
axisValueList{:} を，組み合わせを列挙する combinations に渡しています．splitTableByAxes(T, ["N", "train"]); の場合，結果は次のような table になります．
N train
__ _____
1 true
1 false
2 true
2 false
...
各行が一つの条件に対応しています．二重 for ループで nIdx と trainIdx を動かしていた処理を，条件の一覧表として先に作っていると考えると分かりやすいでしょう．
各条件について filterTable 関数を適用する
次に，axisCombinations を一行ずつ取り出し，変数名とともにセル配列 filterArgs に格納します．その後，これを filterTable に入力して条件を満たす行を抽出します．
% 各条件を満たすテーブルを抽出
Tsplit = cell(height(axisCombinations), 1);
for combinationIdx = 1:height(axisCombinations)
filterArgs = [
num2cell(axisVars)
table2cell(axisCombinations(combinationIdx, :))
];
Tsplit{combinationIdx} = ...
filterTable(T, filterArgs{:});
end
N = 3，train = false に対応する行では，filterArgs は次のようになります．
filterArgs = {
"N", "train"
3, false
};
MATLAB は cell 配列を列優先で展開するため，filterArgs{:} は "N", 3, "train", false に対応します．したがって，以下の二つが一致します．
filterTable(T, filterArgs{:})
filterTable(T, "N", 3, "train", false)
cell 配列の各次元を軸に対応させる
この時点の Tsplit は，組み合わせを縦に並べた cell 配列です．最後に，1次元目が axisVars(1)，2次元目が axisVars(2)，という形になるように形を整えます．
% セル配列の各次元を軸に対応させる
Tsplit = reshape(Tsplit, [fliplr(axisSizes), 1]); % 多次元配列に変形
Tsplit = permute(Tsplit, [numAxes:-1:1, numAxes + 1]); % 軸の順番を逆転
combinations の出力では，後ろの軸ほど速く変化します．一方，MATLAB の reshape は1次元目から順番に要素を詰めます．そのため，まず axisSizes を反転した形に reshape し，その後 permute で次元の順番を戻しています．
reduceTableByAxes の実装
完成版
function [reducedValues, axisValues] = ...
reduceTableByAxes(T, axisVars, reduceFcn, dataVar)
arguments
T table
axisVars (1,:) string
reduceFcn (1,1) function_handle
dataVar (1,1) string
end
[Tsplit, axisValues] = splitTableByAxes(T, axisVars);
% 各グループに関数を適用
fun = @(T) reduceFcn(T.(dataVar));
reducedValues = cellfun(fun, Tsplit, UniformOutput=false);
% 出力が空配列またはスカラーであることを確認
isValid = cellfun(@(x) isempty(x) || isscalar(x), reducedValues);
if ~all(isValid, "all")
error("reduceTableByAxes:InvalidOutputSize", ...
"reduceFcn must return an empty array or a scalar for each group.");
end
% 空の出力を NaN に置き換えて数値配列へ変換
reducedValues(cellfun(@isempty, reducedValues)) = {NaN};
reducedValues = cell2mat(reducedValues);
end
コードの解説はこちら
cellfun で各グループを集約する
splitTableByAxes によって得られる Tsplit はセル配列で，その各要素はテーブルです．したがって，for 文でその各要素を取り出し，その Tij.error の平均値を求めればよいことになります．
for i = 1:10
for j = 1:2
Tij = Tsplit{i,j};
averageError(i,j) = mean(Tij.error);
end
end
cellfun を使うことで，このような「セル配列の各要素に同一の関数を適用する操作」をより簡潔に書くことができます．上記の例の場合，二重の for 文部分を以下の2行に置き換えることができます．
fun = @(T) mean(T.error); % T を受け取って error の平均値を返す無名関数
reducedValues = cellfun(fun, Tsplit);
reduceTableByAxes は，値を集約する関数（ここでは @mean）と，計算対象となる変数（ここでは error）を入力として指定することで汎用性を高めたものです．
空のグループを含む場合の問題
上記の素朴な実装は，fun がすべてのグループに対してスカラーを返す場合には正しく動作します．しかし，splitTableByAxes は各軸に含まれる値の直積をとるため，元の table に存在しない条件の組み合わせに対応する要素は空の table になります．このとき，reduceFcn として @max を指定すると，空のグループに対する出力は空の配列になり，エラーが生じます．
そこで，まず UniformOutput=false を指定し，各グループに対する reduceFcn の出力を cell 配列として受け取ります．この段階では，スカラーと空配列が混在していてもエラーになりません．
fun = @(T) reduceFcn(T.(dataVar));
reducedValues = cellfun(fun, Tsplit, UniformOutput=false);
続いて，各要素が空配列またはスカラーであることを確認します．2要素以上を含む配列が返された場合，集約結果を元の軸構造をもつ一つの配列にまとめられないため，エラーを出します．
isValid = cellfun(@(x) isempty(x) || isscalar(x), reducedValues);
if ~all(isValid, "all")
error("reduceTableByAxes:InvalidOutputSize", ...
"reduceFcn must return an empty array or a scalar for each group.");
end
最後に，空の要素を {NaN} に置き換え，cell2mat で通常の配列へ変換します．これにより，@max のように空入力に対して空配列を返す関数もそのまま指定できます．
reducedValues(cellfun(@isempty, reducedValues)) = {NaN};
reducedValues = cell2mat(reducedValues);
まとめ
本記事では，筆者が普段利用している，table データを整理するための関数を紹介しました．これらの関数を活用することで多重の for 文や論理インデックスの書き直しを削減し，table データの整理を効率化できるはずです！
付録：ダミーデータ作成用関数
function T = dummyData
rng(1)
x = (-2:0.2:2)';
y_observe = x.^3 - 2 * x + 0.3 * randn(size(x));
train = rand(size(x)) < 0.7;
test = ~train;
x_train = x(train);
y_train = y_observe(train);
x_test = x(test);
y_test = y_observe(test);
fit1 = polyfit(x_train, y_train, 1);
fit3 = polyfit(x_train, y_train, 3);
fit10 = polyfit(x_train, y_train, 10);
figure; hold on;
fplot(@(x) polyval(fit1, x), [-2.5 2.5], "k:", "DisplayName", "N=1");
fplot(@(x) polyval(fit3, x), [-2.5 2.5], "k-", "DisplayName", "N=3");
fplot(@(x) polyval(fit10, x), [-2.5 2.5], "k--", "DisplayName", "N=10");
scatter(x_train, y_train, "ro", "DisplayName", "train=true");
scatter(x_test, y_test, "bx", "DisplayName", "train=false");
xlabel("x");
ylabel("y");
legend("Location","northwest");
ylim([-5 5]);
grid on;
T = table();
for n = 1:10
y_predict = polyval(polyfit(x_train, y_train, n), x);
error = abs(y_observe - y_predict);
N = repmat(n, numel(x), 1);
T = [T; table(x, y_observe, y_predict, error, train, N)];
end
hold off
end