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

    function [h,Y_Pred,YRF] = plot(varargin)
%PLOT plot response for two-stage model
%
%  [h,Y_Pred,YRF] = plot(TS,<TS2,TS3,...,TSn>,X,Y,PlotOpts,axhand );

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

IsModel= cellfun('isclass',varargin,'xregtwostage');
AllTSModels= varargin(IsModel);

Args= varargin(~IsModel);

% Data
X= Args{1};

Xs= double(X{1});
XG= X{2};

Ys= Args{2};

if length(Args)<3
    Args{3}= [];
end
[~,TransformPlot,CIFlag,AbsX,ModelRange,CalcPRESS]= iGetPlotOpts(Args{3});

if length(Args)==4 && isgraphics(Args{4})
    AxHand= Args{4};
else
    AxHand= gca;
end



NumModels= length(AllTSModels);

TS= AllTSModels{1};
nl= nlfactors(TS);

doShowRF = nl==1 && size(XG,1)==1;

% find default degrees of freedom
[df1,ni1]= iCILevel(AllTSModels{1},CIFlag);


% set list of styles for validation tool mulitselect
MarkerStyles= iMarkerStyles(AxHand,NumModels);

% cell arrays for building plots
fitPlot= cell(3*length(AllTSModels),1);
pointPlot= fitPlot;
CIPlot= fitPlot;


Y_Pred= cell(1,length(AllTSModels));
YRF= cell(1,length(AllTSModels));


for ModNo= 1:NumModels;
    % loop 
    
    TS= AllTSModels{ModNo};
    
    
    [Ypoints,Yrf,Datum,L]= iEval(TS,X,CalcPRESS);
    
    Y_Pred{ModNo}= Ypoints;
    YRF{ModNo}= Yrf;

    % make the ylabel reflect the
    Trans= TransformPlot && HasTransform(L);
    
    if Trans
        Ydata= ytrans(L,double(Ys));
        Ypoints= ytrans(L,Ypoints);
    else
        Ydata= double(Ys);
    end
    
    ylabstr= ResponseLabel(L,Trans);
    
    IsMultiInput = nl>1 && all(InputFactorTypes(L)==1);
    
    
    [LB,UB]= iCalcRange(L,Xs,Ydata,IsMultiInput,ModelRange);
    
    
    if IsMultiInput
        % plot y=x as solid line
        xline= [LB UB];
        yline= xline;
        
        
        Xpoints= Ydata;
        
        xlabstr= ylabstr;

    else
        Xpoints= Xs(:,1);
        
        xlabstr= InputLabels(TS);
        if nl==1
            % plot solid response line
            xline= linspace(LB,UB,100)';
            yline= L(xline);
            if Trans
                yline= ytrans(L,yline);
            end
        else
            xline= Xpoints;
            yline= Ypoints;
        end

        if ~AbsX
            % take datum away from input 1 to show plot relative to datum
            xline= xline-Datum;
            Xpoints= Xpoints-Datum;
            xlabstr = [xlabstr{1} ' - DATUM'];
        else
            % raw x data
            xlabstr= xlabstr{1};
        end
    end
    
    % main line
    fitPlot(3*ModNo-2:3*ModNo)= {xline , yline , '-'};
    
    % Prediction Points
    pointPlot(3*ModNo-2:3*ModNo)= {Xpoints , Ypoints , MarkerStyles{ModNo} };
    
    % Confidence Interval Plots
    if CIFlag && pevcheck(TS)
        [~,ni]= iCILevel(TS,CIFlag,df1,ni1);
        
        [xci,yci]= iConfidenceInterval(TS,{Xs,XG},Xpoints,Ypoints,IsMultiInput,Trans,ni,UB-LB);
        CIPlot(3*ModNo-2:3*ModNo)= {xci(:) , yci(:) , '-'};
    else
        CIPlot(3*ModNo-2:3*ModNo)= {NaN , NaN , '-'};
    end
    
    
    if doShowRF
        iPlotRF(L,AxHand,Trans,AbsX);
    end
end

% draw plots
h= iDoPlots(fitPlot,pointPlot,CIPlot,AxHand,CIFlag,MarkerStyles);

set(get(AxHand,'Title'),'String',sprintf('Test %3g',testnum(X{1})),...
    'FontWeight','Bold','FontSize',8)
set(get(AxHand,'XLabel'),...
    'String',xlabstr,...
    'Interpreter','none',...
    'FontSize',8);
set(get(AxHand,'YLabel'),...
      'String',ylabstr,...
      'Interpreter','none',...
      'FontSize',8);



%--------------------------------------------------------------------------
function hPoints= iDoPlots(fitPlot,pointPlot,CIPlot,AxHand,CIFlag,MarkerStyles)

% calls MATLAB plot with fitplot={xdata,ydata,linestyle}

ColorIndex = AxHand.ColorOrderIndex;
h=plot(fitPlot{:},'Parent',AxHand);
set(h,'Tag','twostagefit');
xregGui.setLegendData(h, 'Model fit');
% make sure next plots have the same colors as the first
AxHand.ColorOrderIndex = ColorIndex;
if length(MarkerStyles)==1
    hPoints= plot(pointPlot{:},'Parent',AxHand,'MarkerSize',15);
else
    hPoints= plot(pointPlot{:},'Parent',AxHand,'MarkerSize',5);
end
set(hPoints,'Tag','twostagepoints');
xregGui.setLegendData(hPoints, 'Model predictions');

if CIFlag
    % make sure next plots have the same colors as the first
    AxHand.ColorOrderIndex = ColorIndex;
    hCI=plot(CIPlot{:},'Parent',AxHand,'LineWidth',1);
    set(hCI,'Tag','twostageCI');
    xregGui.setLegendData(hCI, false);
end


%--------------------------------------------------------------------------
function [bdflag,Trans,CIFlag,AbsX,ModelRange,CalcPRESS]= iGetPlotOpts(PlotOpts)

% plotopts must be a row vector
if (islogical(PlotOpts) || isnumeric(PlotOpts)) && size(PlotOpts,1)==1 && (numel(PlotOpts)==5 || numel(PlotOpts)==6)
    PlotOpts= num2cell(PlotOpts);
    if numel(PlotOpts)==5
        % CalcPRESS
        PlotOpts{6}= 0;
    end
else
    PlotOpts   = {0,0,0,1,0,0};
end

[bdflag,Trans,CIFlag,AbsX,ModelRange,CalcPRESS]= deal(PlotOpts{:});
if CalcPRESS
    % turn off CI's
    CIFlag=0;
end

%--------------------------------------------------------------------------
function [xci,yci]= iConfidenceInterval(TS,Xdata,Xpoints,Ypoints,IsMultiInput,Trans,ni,Range)

xci= Xpoints;
XLci= code(TS,double(Xdata{1}),1:size(Xdata{1},2));
x2= gcode(TS,double(Xdata{2}));
% use transform
if size(x2,1)==size(XLci,1)
    ts= ni*sqrt(pev(TS,[XLci,x2],0,~Trans));
elseif IsMultiInput
    ts= ni*sqrt(pev(TS,{XLci,x2},0,~Trans));
else
    ts= ni*sqrt(pevgrid(TS,[num2cell(XLci,1),num2cell(x2)],0,~Trans));
end

N= size(xci,1);

xci=xci(:,1);
% make confidence interval lines
nanM= NaN(N,1);
if N>=100
    % single bar for large samples
    yci= Ypoints(:,ones(3,1))'+[-ts ts nanM]';
    if IsMultiInput
        [xci,ind] = sort(xci);
        xci= [xci xci nanM]';
        yci= yci(:,ind);
    else
        xci= [xci xci nanM]';
    end	
else
    % bars with tops and tails if less than 100
    dx= Range/150;
    yci= Ypoints(:,ones(9,1))'+[-ts ts nanM ts ts nanM -ts -ts nanM]';
    if IsMultiInput
        [xci,ind] = sort(xci);
        xci= [xci xci nanM xci-dx xci+dx nanM xci-dx xci+dx nanM]';
        yci= yci(:,ind);
    else
        xci= [xci xci nanM xci-dx xci+dx nanM xci-dx xci+dx nanM]';
    end	
end



%--------------------------------------------------------------------------
function [LB,UB]= iCalcRange(L,Xdata,Ydata,IsMultiInput,ModelRange)

if IsMultiInput
    LB= min(Ydata);
    UB= max(Ydata);
else
    % calculate bounds for plots from rf values and Xs(:,1)
    fVals= unique(get(L,'Values'));

    if ModelRange
        Bnds= getcode(L);
        if ~((Bnds(1,1)==-1 || Bnds(1,1)==0)  && Bnds(1,2)==1) 
            LB= Bnds(1,1);
            UB= Bnds(1,2);
        else
            % Old Model with no ranges
            LB= min(Xdata(:,1));
            UB= max(Xdata(:,1));
        end   
    elseif any(fVals)
        if DatumType(L)
            D = datum(L);
        else
            D = 0;
        end
        LB= min([fVals(:,1)+D(1);Xdata(:,1)]);
        UB= max([fVals(:,1)+D(1);Xdata(:,1)]);
    else
        LB= min(Xdata(:,1));
        UB= max(Xdata(:,1));
    end
end


%--------------------------------------------------------------------------
function iPlotRF(L,AxHand,Trans,AbsX)
        
if DatumType(L)
    D = datum(L);
else
    D = 0;
end
fVals= unique(get(L,'Values'));
featvals= L(fVals+D);
datval= L(D);

if Trans
    featvals= ytrans(L,featvals);
    datval= ytrans(L,datval);
end

if ~AbsX
    fpts= fVals;
    dpt= 0;
else
    dpt= D;
    fpts= fVals+D;
end
% plot response feature values 
hL = line('XData',fpts(fVals~=0),'YData',featvals(fVals~=0),...
    'Marker','+','Color','m','LineStyle','none',...
    'Tag','twostagefcnval','Parent',AxHand);
xregGui.setLegendData(hL, 'Response features');
if DatumType(L)
    % plot datum if required
    hL = line('XData',dpt,'YData',datval,...
        'Marker','*','Color','m','LineStyle','none',...
        'Tag','twostagedatum','Parent',AxHand);
    xregGui.setLegendData(hL, 'Datum');
end


%--------------------------------------------------------------------------
function [Y,Yrf,Datum,L]= iEval(TS,Data,CalcPRESS)

if CalcPRESS
    [Y,Yrf,Datum,p]= presspred(TS,Data,CalcPRESS);
else
    [Y,Yrf,Datum,p]= EvalModel(TS,Data);
end

if size(Yrf,1)>1
    Yrf= mean(Yrf,1);
    Datum= mean(Datum,1);
    p= mean(p,1);
end

% form local model
L= get(TS,'local');
if DatumType(L)
    L= datum(L,Datum);
end
L= update(L,p,[]);


%--------------------------------------------------------------------------
function [df,ni]= iCILevel(TS,CIFlag,df1,ni1)

df=dferror(TS);
if nargin>2 && df1==df
    ni= ni1;
else
    
    if CIFlag==1
        alpha= 0.975;
    else
        alpha= 1-(1-CIFlag/100)/2;
    end
    
    % find default degrees of freedom
    if ~isfinite(df) || df > 200
        df= Inf;
        ni = norminv(alpha);
    else
        ni = tinv(alpha,df);
    end
end

%--------------------------------------------------------------------------
function MarkerStyles= iMarkerStyles(AxHand,NumModels)

% set list of styles for validation tool mulitselect
MarkerStyles= get(AxHand,'LineStyleOrder');
% if plotting regular twostage at local node, want dots, not ML default '-'
if ~iscell(MarkerStyles)
    MarkerStyles= {'.'};
end

while length(MarkerStyles)<NumModels;
    MarkerStyles=[MarkerStyles;MarkerStyles];
end