%% Plot Confidence Intervals of PK Model Predictions
% 
%% Load Data
% Load the sample data to fit.
load data10_32R.mat
gData = groupedData(data);
gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter'};
sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},'Marker','+',...
            'LineStyle','none');
%% Create Model
% Create a two-compartment model.
pkmd                 = PKModelDesign;
pkc1                 = addCompartment(pkmd,'Central');
pkc1.DosingType      = 'Infusion';
pkc1.EliminationType = 'linear-clearance';
pkc1.HasResponseVariable = true;
pkc2                 = addCompartment(pkmd,'Peripheral');
model                = construct(pkmd);
configset            = getconfigset(model);
configset.CompileOptions.UnitConversion = true;
%% Define Dosing
% Define the infusion dose. 
dose             = sbiodose('dose','TargetName','Drug_Central');
dose.StartTime   = 0;
dose.Amount      = 100;
dose.Rate        = 50;
dose.AmountUnits = 'milligram';
dose.TimeUnits   = 'hour';
dose.RateUnits   = 'milligram/hour';
%% Define Parameters
% Define parameters to estimate. 
responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};
paramsToEstimate   = {'log(Central)','log(Peripheral)','Q12','Cl_Central'};
estimatedParam     = estimatedInfo(paramsToEstimate,...
                                   'InitialValue',[1 1 1 1],...
                                   'Bounds',[0.1 3;0.1 10;0 10;0.1 2]);
%% Fit Model
% Perform an unpooled fit, that is, one set of estimated parameters for
% each patient.
unpooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',false);
%% Compute Confidence Intervals for Model Predictions
% Compute 95% confidence intervals for predicted model responses in the
% unpooled fit using the Gaussian approximation.
ciPredUnpooled = sbiopredictionci(unpooledFit);
%% 
% Plot the confidence intervals. If the estimation status of a confidence
% interval is |constrained| or |not estimable|, the function uses the second
% default color (red). Otherwise, the function uses the first default color
% (blue). To see the color order, type
% |get(groot,'defaultAxesColorOrder')|.
%%
% Each group is displayed in each column, from left to right, in the same
% order that they appear in the |GroupNames| property of the object, which
% is used to label each column.
plot(ciPredUnpooled)
%%
% Plot using a single color (black).
plot(ciPredUnpooled,'Color',[0 0 0])

%% 
% Copyright 2012 The MathWorks, Inc.