www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregmodel/pevgrid.m

    function [PEV,Xout,Xg,Y] = pevgrid(m,Values,Natural,varargin)
%PEVGRID Evaluate Prediction Error Variance for model over grid
%
%  PEV = evalpev(m,Values,Natural)
%    m is the model. InitStore must be called on m before this function
%    Values cell array defining grid e.g. Values = {-1:.1:1,0,-1:.1:1}
%    Natural==1 if Values is in natural units
%
%  If y data is available, PEV = x'* s*inv(X'*X) * x, otherwise PEV = x'*
%  inv(X'*X) * x.  NDGRID is used to define the grid. The shape of PEV is
%  the same as the grid shape returned by NDGRID.
%
%  See also XREGLINEAR/EVALPEV, NDGRID

%  Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc.


if nargin<=2
    Natural=1;
end
xc= Values;
s= cellfun('prodofsize',xc);
MaxSize= min(10000,prod(s));
dimGrid = sum(s>1);

if Natural
    for i=1:length(xc);
        % do the coding on each entry individually as this is
        % computationally much cheaper
        xc{i}= code(m,xc{i}(:),i);
    end
end
n= length(xc);
IT2=  ~all(InputFactorTypes(m)==1);
if n>1

    X=Values;
    if IT2
        Xout= Values;
        % transient data dependent with time dependent outputs
        X= xc;
        MaxSize= length(X{1});
    else
        if dimGrid>1

            % Generate N-D grid for evaluation
            [X{s>1}] = ndgrid(xc{s>1});

            if nargout>1
                [Xout{s>1}] = ndgrid(Values{s>1});
            end
        else
            X= xc;
            Xout= Values;
        end
    end
else
    X= xc;
    Xout= Values;
end





Xg =zeros(MaxSize,n);
for i= find(s==1)
    % setup scalar inputs
    Xg(:,i)= xc{i};
end


if IT2
    Nevals= 1;
    PEV= zeros(numel(X{1}),1);
    if nargout>3
        Y= PEV;
    end
else
    Y = zeros(prod(s),1);
    Nevals = floor(prod(s)/MaxSize);
    PEV = zeros(prod(s),1);
    if nargout>3
        Y = PEV;
    end
end

for i= 1:Nevals
    % do evaluations in blocks of MaxSize
    ind= (MaxSize*(i-1)+1:MaxSize*i)';
    for j=find(s>1)
        % non scalar inputs
        Xg(:,j)= X{j}(ind);
    end
    if nargout>3
        [PEV(ind),Y(ind)]= pev(m,Xg,0,varargin{:});
    else
        PEV(ind)= pev(m,Xg,0,varargin{:});
    end
end
if ~IT2 && MaxSize*i<prod(s)
    % last block
    ind= (floor(prod(s)/MaxSize)*MaxSize+1:prod(s))';
    Xg= Xg(1:length(ind),:);
    for j=find(s>1)
        % non scalar inputs
        Xg(:,j)= X{j}(ind);
    end
    if nargout>3
        [PEV(ind),Y(ind)]= pev(m,Xg,0,varargin{:});
    else
        PEV(ind)= pev(m,Xg,0,varargin{:});
    end
end

% reshape table
if length(s)>1 && ~IT2
    PEV = reshape(PEV,s);
    if nargout>3
        Y = reshape(Y,s);
    end
end

if nargout>2
    [dum,d] = max(s);
    Xg = zeros(numel(X{d}),n);
    for i=1:n
        Xg(:,i) = X{i}(:);
    end
end