MapReduce を使用するロジスティック回帰モデルの近似

ライブスクリプトを開く

この例では、mapreduce を使用して、単一の予測子を使用して簡単なロジスティック回帰を実行する方法を示します。複数の mapreduce 呼び出しを連結して反復アルゴリズムを実行する方法も示します。各反復では独立したデータパスが必要なため、無名関数は 1 つの反復から次の反復に情報を渡して、マッパーに直接情報を提供します。

データの準備

airlinesmall.csv データセットを使用してデータストアを作成します。この 12 MB のデータセットには、到着時間と出発時間を含む、いくつかの航空会社のフライト情報が 29 列に含まれます。この例では、対象の変数は ArrDelay (フライトの到着遅延時間) と Distance (総飛行距離) です。

ds = tabularTextDatastore('airlinesmall.csv', 'TreatAsMissing', 'NA');
ds.SelectedVariableNames = {'ArrDelay', 'Distance'};

データストアは、既定では 'NA' 値を欠損として扱い、欠損値を NaN 値に置換します。さらに、SelectedVariableNames プロパティを使用すると、指定された対象の変数のみを処理して、preview を使用して確認できます。

preview(ds)

ans=8×2 table
    ArrDelay    Distance
    ________    ________

        8         308   
        8         296   
       21         480   
       13         296   
        4         373   
       59         308   
        3         447   
       11         954

ロジスティック回帰の実行

ロジスティック回帰は、あるイベントの確率を他の変数の関数としてモデル化する方法です。この例では、ロジスティック回帰により、フライトが 20 分より遅延する確率を、1,000 マイル単位の飛行距離の関数としてモデル化します。

このロジスティック回帰を実行するために、map 関数と reduce 関数は集合的に、現在の係数の値に基づいて重み付き最小二乗回帰を実行します。マッパーは、入力データのブロックごとに重み付き二乗和とクロス積を計算します。

map 関数のファイルを表示します。

function logitMapper(b,t,~,intermKVStore)
  % Get data input table and remove any rows with missing values
  y = t.ArrDelay;
  x = t.Distance;
  t = ~isnan(x) & ~isnan(y);
  y = y(t)>20;                 % late by more than 20 min
  x = x(t)/1000;               % distance in thousands of miles

  % Compute the linear combination of the predictors, and the estimated mean
  % probabilities, based on the coefficients from the previous iteration
  if ~isempty(b)
    % Compute xb as the linear combination using the current coefficient
    % values, and derive mean probabilities mu from them
    xb = b(1)+b(2)*x;
    mu = 1./(1+exp(-xb));
  else
    % This is the first iteration. Compute starting values for mu that are
    % 1/4 if y=0 and 3/4 if y=1. Derive xb values from them.
    mu = (y+.5)/2;
    xb = log(mu./(1-mu)); 
  end

  % To perform weighted least squares, compute a sum of squares and cross
  % products matrix:
  %      (X'*W*X) = (X1'*W1*X1) + (X2'*W2*X2) + ... + (Xn'*Wn*Xn),
  % where X = [X1;X2;...;Xn]  and  W = [W1;W2;...;Wn].
  %
  % The mapper receives one chunk at a time and computes one of the terms on
  % the right hand side. The reducer adds all of the terms to get the
  % quantity on the left hand side, and then performs the regression.
  w = (mu.*(1-mu));                  % weights
  z = xb + (y - mu) .* 1./w;         % adjusted response

  X = [ones(size(x)),x,z];           % matrix of unweighted data
  wss = X' * bsxfun(@times,w,X);     % weighted cross-products X1'*W1*X1

  % Store the results for this part of the data.
  add(intermKVStore, 'key', wss);
end

リデューサーは、二乗和とクロス積から回帰係数の推定値を計算します。

reduce 関数のファイルを表示します。

function logitReducer(~,intermValIter,outKVStore)
  % We will operate over chunks of the data, updating the count, mean, and
  % covariance each time we add a new chunk
  old = 0;

  % We want to perform weighted least squares. We do this by computing a sum
  % of squares and cross products matrix
  %      M = (X'*W*X) = (X1'*W1*X1) + (X2'*W2*X2) + ... + (Xn'*Wn*Xn)
  % where X = X1;X2;...;Xn]  and  W = [W1;W2;...;Wn].
  %
  % The mapper has computed the terms on the right hand side. Here in the
  % reducer we just add them up.

  while hasnext(intermValIter)
    new = getnext(intermValIter);
    old = old+new;
  end
  M = old;  % the value on the left hand side

  % Compute coefficients estimates from M. M is a matrix of sums of squares
  % and cross products for [X Y] where X is the design matrix including a
  % constant term and Y is the adjusted response for this iteration. In other
  % words, Y has been included as an additional column of X. First we
  % separate them by extracting the X'*W*X part and the X'*W*Y part.
  XtWX = M(1:end-1,1:end-1);
  XtWY = M(1:end-1,end);

  % Solve the normal equations.
  b = XtWX\XtWY;

  % Return the vector of coefficient estimates.
  add(outKVStore, 'key', b);
end

mapreduce の実行

mapreduce の呼び出しをループ内に囲むことによって、mapreduce を繰り返し実行します。ループは、反復回数 5 回を最高として、収束条件が満たされるまで実行されます。

% Define the coefficient vector, starting as empty for the first iteration.
b = []; 

for iteration = 1:5
    b_old = b;
    iteration
    
    % Here we will use an anonymous function as our mapper. This function
    % definition includes the value of b computed in the previous
    % iteration.
    mapper = @(t,ignore,intermKVStore) logitMapper(b,t,ignore,intermKVStore);
    result = mapreduce(ds, mapper, @logitReducer, 'Display', 'off');
    
    tbl = readall(result);
    b = tbl.Value{1}
    
    % Stop iterating if we have converged.
    if ~isempty(b_old) && ...
       ~any(abs(b-b_old) > 1e-6 * abs(b_old))
       break
    end
end

iteration = 1

iteration = 2

iteration = 3

iteration = 4

結果の表示

結果の回帰係数の推定値を使用して、確率曲線をプロットします。この曲線は、フライトの遅延時間が 20 分を超える確率を、飛行距離の関数として示します。

xx = linspace(0,4000);
yy = 1./(1+exp(-b(1)-b(2)*(xx/1000)));
plot(xx,yy); 
xlabel('Distance');
ylabel('Prob[Delay>20]')

Figure contains an axes object. The axes object contains an object of type line.

ローカル関数

ここに挙げるのは、mapreduce がデータに適用する map 関数と reduce 関数です。

function logitMapper(b,t,~,intermKVStore)
  % Get data input table and remove any rows with missing values
  y = t.ArrDelay;
  x = t.Distance;
  t = ~isnan(x) & ~isnan(y);
  y = y(t)>20;                 % late by more than 20 min
  x = x(t)/1000;               % distance in thousands of miles

  % Compute the linear combination of the predictors, and the estimated mean
  % probabilities, based on the coefficients from the previous iteration
  if ~isempty(b)
    % Compute xb as the linear combination using the current coefficient
    % values, and derive mean probabilities mu from them
    xb = b(1)+b(2)*x;
    mu = 1./(1+exp(-xb));
  else
    % This is the first iteration. Compute starting values for mu that are
    % 1/4 if y=0 and 3/4 if y=1. Derive xb values from them.
    mu = (y+.5)/2;
    xb = log(mu./(1-mu)); 
  end

  % To perform weighted least squares, compute a sum of squares and cross
  % products matrix:
  %      (X'*W*X) = (X1'*W1*X1) + (X2'*W2*X2) + ... + (Xn'*Wn*Xn),
  % where X = [X1;X2;...;Xn]  and  W = [W1;W2;...;Wn].
  %
  % The mapper receives one chunk at a time and computes one of the terms on
  % the right hand side. The reducer adds all of the terms to get the
  % quantity on the left hand side, and then performs the regression.
  w = (mu.*(1-mu));                  % weights
  z = xb + (y - mu) .* 1./w;         % adjusted response

  X = [ones(size(x)),x,z];           % matrix of unweighted data
  wss = X' * bsxfun(@times,w,X);     % weighted cross-products X1'*W1*X1

  % Store the results for this part of the data.
  add(intermKVStore, 'key', wss);
end
%-----------------------------------------------------------------------------
function logitReducer(~,intermValIter,outKVStore)
  % We will operate over chunks of the data, updating the count, mean, and
  % covariance each time we add a new chunk
  old = 0;

  % We want to perform weighted least squares. We do this by computing a sum
  % of squares and cross products matrix
  %      M = (X'*W*X) = (X1'*W1*X1) + (X2'*W2*X2) + ... + (Xn'*Wn*Xn)
  % where X = X1;X2;...;Xn]  and  W = [W1;W2;...;Wn].
  %
  % The mapper has computed the terms on the right hand side. Here in the
  % reducer we just add them up.

  while hasnext(intermValIter)
    new = getnext(intermValIter);
    old = old+new;
  end
  M = old;  % the value on the left hand side

  % Compute coefficients estimates from M. M is a matrix of sums of squares
  % and cross products for [X Y] where X is the design matrix including a
  % constant term and Y is the adjusted response for this iteration. In other
  % words, Y has been included as an additional column of X. First we
  % separate them by extracting the X'*W*X part and the X'*W*Y part.
  XtWX = M(1:end-1,1:end-1);
  XtWY = M(1:end-1,end);

  % Solve the normal equations.
  b = XtWX\XtWY;

  % Return the vector of coefficient estimates.
  add(outKVStore, 'key', b);
end
%-----------------------------------------------------------------------------

参考

mapreduce | tabularTextDatastore