How to organize plot of sbiosobol results
3 次查看(过去 30 天)
显示 更早的评论
Hello,
I used to use the GlobalSensitivityAnalysisApp from file exhange. However, I would like to use the script to carry out the sbiosobol in order to be able to define more details for my global sensitivity analysis (sobolResults = sbiosobol(...)). At the very bottom of the GlobalSensitivityAnalysisApp user interface, there is a plot setting for defining how many parameters and obserables per plot. How could I do the same plot setting with code for the result from the sbiosobol?
Please give me any suggestion or advice. Thank you very much.
Have a nice day.
Jesse
0 个评论
采纳的回答
Florian Augustin
2021-11-10
Hello Jesse,
You can use the name-value arguments "Parameters" and "Observables" of the plot and bar functions to specify which sensitivity inputs (parameters) and sensitivity outputs (observables) are plotted. You can specify the values as names of input parameters / observables, or as numeric indices. The plot and bar function respect the order of the specified parameters / observables / indices, this way you can also sort the Sobol indices.
For example, assuming you have 6 sensitivity inputs / parameters and want to split the plot into two figures, you can run:
% sobolResults = sbiosobol(...);
plot(sobolResults, "Parameters", 1:3);
plot(sobolResults, "Parameters", 4:6);
Best,
Florian
6 个评论
Florian Augustin
2021-11-15
Hi Jesse,
You should be able to plot on top of any GSA plot by turning on "hold" for the axes you want to plot into. You'd first have to get the axes from the plot, hold the existing data of the axes, and then plot on top of it:
% sobolResults = sbiosobol(...);
hFig = plotData(sobolResults); % get figure handle
hAxs = findobj(hFig, "Type", "Axes") % find axes to draw on
selectedAxes = hAxs(1); % assuming axes 1 is the one you want to add content to
hold(selectedAxes, "on"); % hold existing content in axes
plot(selectedAxes, ...); % plot into axes
hold(selectedAxes, "off");
The function sbiompgsa only evaluates the whole classifier and does not store individual evaluations of the left- and right-hand sides. Understanding what is going on during the evaluation of classifiers can be tricky; let me share my workflow how I ususally "debug" the evaluation of classifiers. I typically run a couple of model simulations to see the general ballpark of the ranges of the left- and right-hand sides of the classifier's inequality.
% sobolResults = sbiosobol(...);
% Get SimFunction to reproduce model simulations
simFun = sobolResults.SimulationInfo.SimFunction;
% Get time steps
simTimes = sobolResults.Time;
% Get samples for parameter inputs
paramSamples = sobolResults.ParameterSamples{:,:};
% This is equivalent to writing:
% paramSamples = getSimulationResults(sobolResults, 1);
% If simulating all paramSamples takes too long, you can just select a subset, e.g.
% every 10th sample:
% paramSamples = sobolResults.ParameterSamples{1:10:end, :};
% Run model simulations:
simData = simFun(paramSamples, [], [], simTimes);
% Compute left- and right-hand side of inequality of classifier:
simData = addobservable(simData, ["lhs", "rhs"], ...
["abs(trapz(time, PlasmaConc) - AUC_plasma_mean)", ...
"0.5*AUC_plasma_mean"]);
classifierData = selectbyname(simData, ["lhs", "rhs"]);
% Plot lhs and rhs
sbioplot(classifierData);
SimBiology's SimFunctions and how to use them is documented here. The above sample code should give you a spaghetti plot of values of the left- and right-hand side of the classifiers for a range of parameter samples. I assume in your case AUC_plasma_mean is a constant scalar parameter, so you should be see from the evaluation of all lhs evaluations are above, below, or distributed around the rhs. This should match the evaluation of the classifier, i.e. the property SupportHypothesis, on the results object returned from sbiompgsa.
No worries about having a public conversation on this tread. You raise very interesting questions other SimBiology users may also be interested in. However, feel free to send me an email via my community profile if you have model specific questions. I know many users are unable to share model details publicly.
I hope this helps.
Best,
Florian
更多回答(1 个)
Jesse Chao
2021-11-18
1 个评论
Florian Augustin
2021-11-19
Hi Jesse,
The SimFunction you get from sobolResults already knows about the model and variants you specified in the call to sbiosobol. But I think the problem is that the dosing information is missing. My example assumed that there are no doses. To apply doses, you need to supply them as follows:
% Assuming doses is a row vector of SimBiology RepeatDose / ScheduleDose objects.
% If you are not sure, you can turn your doses into a row vector as follows:
% doses = doses(:)';
simData = simFun(paramSamples, [], getTable(doses), simTimes);
The reason you need to supply the doses is because they describe "runtime" information that you can change between model simulations. All the different ways to apply doses in SimFunctions is documented here (look for the input argument "u").
The first step would be to get the model simulation working, then we can look at the observables for the left- and right-hand sides of the classifier. Feel free to contact me via my profile page.
I hope that helped.
-Florian
社区
更多回答在 SimBiology Community
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Perform Sensitivity Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!