forwardselect.m

人工神经网络的源码编程

      ldiagHoTHo2 = ldiagHoTHo  .* ldiagHoTHo;
      ldiagHoTHo3 = ldiagHoTHo2 .* ldiagHoTHo;
      YP2Y = YY - sum(l2diagHoTHo .* YTHoTHoY ./ ldiagHoTHo2);
      WAW = sum(YTHoTHoY ./ ldiagHoTHo3);
      tA = sum(1 ./ ldiagHoTHo);
      tA2 = sum(1 ./ ldiagHoTHo2);
      g = sum(diagHoTHo ./ ldiagHoTHo);

    else

      A = inv(Hm' * Hm + lambda * eye(m));
      HY = Hm' * Y;
      W = A * HY;
      PY = Y - Hm * W;
      YP2Y = traceProduct(PY', PY);
      WAW = traceProduct(W', A * W);
      tA = trace(A);
      tA2 = trace(A^2);
      g = m - lambda * tA;

    end

    % exercise different re-estimation methods
    if ReEst == 'u'
      psi = 1 / (2 * (p - g));
    elseif ReEst == 'f'
      psi = p / ((p - g) * (p + g));
    elseif ReEst == 'g'
      psi = 1 / (p - g);
    else
      psi = p * log(p) / (2 * (p - g) * (p + (log(p) - 1) * g));
    end

    % do the re-estimation
    lambda = psi * YP2Y * (tA - lambda * tA2) / WAW;

    % print new value
    if Verbose
      fprintf('%8.3e ', lambda)
    end

    % keep track of multiple lambdas
    if m == 1
      lambdas = lambda;
    else
      lambdas = [lambdas lambda];
    end

  end

  % calculate current score
  if Term == 'o'     % LOO

    if OLS

      % need diag(Pm)
      ldiagHoTHo = lambda + diagHoTHo;
      diagPm = ones(p,1) - diagProduct(Ho, dupCol(1 ./ ldiagHoTHo, p) .* Ho');

      % watch out for zero entries along the diagonal of P
      tooSmall = find(diagPm < eps);
      diagPm(tooSmall) = eps * ones(length(tooSmall),1);

      % need Pm * Y
      PmY = Y - Ho * ((Ho' * Y) ./ dupCol(ldiagHoTHo, n));

      % need inv(diag(diag(Pm))) * Pm * Y
      invDiagPmPmY = PmY ./ dupCol(diagPm, n);

      % compute LOO (Sherman-Morrison-Woodbury formula)
      score = traceProduct(invDiagPmPmY', invDiagPmPmY) / p;

    else

      A = inv(Hm' * Hm + lambda * eye(m));
      HY = Hm' * Y;
      AH = A * Hm';
      PY = Y - Hm * (AH * Y);
      dP = ones(p,1) - diagProduct(Hm, AH);
      dPPY = PY ./ dupCol(dP, n);
      score = traceProduct(dPPY', dPPY) / p;

    end

  else

    % get trace(Pm) and also AH if not using OLS
    if lambda == 0

      if OLS

        tracePm = p - m;

      else

        A = inv(Hm' * Hm);
        AH = A * Hm';
        tracePm = p - m;

      end

    else

      if OLS

        tracePm = p - sum(diagHoTHo ./ (lambda + diagHoTHo));

      else

        A = inv(Hm' * Hm + lambda * eye(m));
        AH = A * Hm';
        tracePm = p - traceProduct(Hm, AH);

      end

    end

    % Get number of effective parameters (g).
    % Watch out for zero trace (usually p = m and lambda = 0).
    if tracePm < eps
      tracePm = eps;
      g = p;
    else
      g = p - tracePm;
    end

    % get mean square error
    if OLS

      if lambda == 0
        YP2Yp = (YY - sum(diagProduct(HoTY,HoTY') ./ diagHoTHo)) / p;
        YP2Yp = (YY - traceProduct(HoTY', HoTY ./ dupCol(diagHoTHo, n))) / p;
      else
        ldiagHoTHo  = lambda + diagHoTHo;
        ldiagHoTHo2 = ldiagHoTHo .^2;
        YP2Yp = (YY - traceProduct(HoTY', HoTY ./ dupCol(ldiagHoTHo, n)) ...
          - lambda * traceProduct(HoTY', HoTY ./ dupCol(ldiagHoTHo2, n))) / p;
      end

    else

      PY = Y - Hm * (AH * Y);
      YP2Yp = traceProduct(PY', PY) / p;

    end

    % get different factors for each method
    if Term == 'e'
      % unexplained variance
      psi = 1;
    elseif Term == 'u'
      % UEV
      psi = p / tracePm;
    elseif Term == 'f'
      % FPE
      psi = (p + g) / tracePm;
    elseif Term == 'g'
      % GCV
      psi = p^2 / tracePm^2;
    else
      % BIC
      psi = (p + (log(p) - 1) * g) / tracePm;
    end

    % finally compute score
    score = psi * YP2Yp;

  end
  if Verbose
    fprintf('%9.3e ', score)
  end

  % are we ready to terminate yet
  if tot >= 1

    Continue = 0;
    if Flops
      fprintf('%8d  ', flops)
    end
    if Verbose
      fprintf('\nvariance all explained ')
    end

  elseif m >= M

    Continue = 0;
    if Flops
      fprintf('%8d  ', flops)
    end
    if Verbose
      fprintf('\nregressors used up ')
    end

  elseif MaxReg > 0 & m >= MaxReg

    Continue = 0;
    if Flops
      fprintf('%8d  ', flops)
    end
    if Verbose
      fprintf('\nlimit of regressors reached ')
    end

  else

    % decide if termination conditions have been met
    if m == 1

      % don't stop here unless threshold being used
      if Term == 'e'

        if (1 - p * score / YY) > Threshold

          Continue = 0;
          if Flops
            fprintf('%8d  ', flops)
          end
          if Verbose
            fprintf('\nexplained variance threshold exceeded ')
          end

        end

      else

        MinScore = score;
        AgeOfMin = 1;

      end       

    else

      % decide between threshold and minimum methods
      if Term == 'e'

        if (1 - p * score / YY) > Threshold

          Continue = 0;
          if Flops
            fprintf('%8d  ', flops)
          end
          if Verbose
            fprintf('\nexplained variance threshold exceeded ')
          end

        end

      else

        % compare old and new score and age of minimum
        if score < MinScore

          % new minimum - don't stop here
          MinScore = score;
          AgeOfMin = 1;

        else

          if AgeOfMin >= MaxAge

            % old minimum has gone on long enough - stop here
            Continue = 0;
            if Flops
              fprintf('%8d  ', flops)
            end
            if Verbose
              fprintf('\nminimum passed ')
            end

	  else

            % old minimum just ages by 1
            AgeOfMin = AgeOfMin + 1;

          end

        end

      end

    end

  end

  if Flops & Continue
    fprintf('%8d  ', flops)
  end
  if Verbose
    fprintf('\n')
  end

end

% don't include last few regressors which aged the minimum
m = m - AgeOfMin;
subset = subset(1:m);

% truncate OLS structures (in case minimum aged and m shrank)
if OLS

  Um = Um(1:m, 1:m);
  diagHoTHo = diagHoTHo(1:m);
  HoTY = HoTY(1:m,:);

end

if nargin > 1

  % actual design matrix being used
  H = F(:, subset);

  if nargout > 2

    % regularisation constant
    if ReEst == 'n'
      l = lambda;
    else
      l = lambdas(m);
    end

    if nargout > 3

      % upper triangular matrix
      if OLS
        U = Um;
      else
        U = eye(m);
      end

      if nargout > 4

        % covariance
        if OLS
          invU = inv(U);
          Ao = diag(1 ./ (diagHoTHo + lambda));
          A = invU * Ao * invU';
        else
          A = inv(H' * H + lambda * eye(m));
        end

        if nargout > 5

          % weight vector (or matrix)
          if OLS
            wo = Ao * HoTY;
            w = invU * wo;
          else
            w = A * (H' * Y);
          end

          if nargout > 6

            % projection matrix
            if OLS
              P = eye(p);
              for j = 1:m
                P = P - Ho(:,j) * Ho(:,j)' / (lambda + diagHoTHo(j));
              end
            else
              P = eye(p) - H * A * H';
            end

          end

        end

      end

    end

  end

end

if Verbose
  fprintf('\n')
end