Formatting for 'addobservables'

3 views (last 30 days)
Shaun Tan
Shaun Tan on 11 Nov 2021
Commented: Shaun Tan on 19 Nov 2021
I've built a SIR epidemic model on SimBiology and want to calculate its sobol indices.
I have 2 parameters, alpha and gamma, which describe infection rate and recovery rate, and 3 species, 'S', 'I' and 'R' to represent the population compartments.
I want to determine peak infections after the simulation (what is the highest value for the 'I' species) and set that as the observable which I will the input into the sbiosobol function to determine effect of alpha and gamma on the value of peak infections.
This is the code:
m1=sbiomodel('SIR');
N=10000;
I0=10;
S=addspecies(m1,'S','InitialAmount',N-I0);
I=addspecies(m1,'I','InitialAmount',I0);
R=addspecies(m1,'R','InitialAmount',0);
r1=addrule(m1,'S = -beta*S*I/N','RuleType','rate');
r2=addrule(m1,'I = beta*S*I/N-gamma*I','RuleType','rate');
r3=addrule(m1,'R = gamma*I','RuleType','rate');
p1=addparameter(m1,'beta','Value',0.2);
p2=addparameter(m1,'gamma','Value',0.1);
p3=addparameter(m1,'N','Value',N);
T=200;
conf=getconfigset(m1,'active');
set(conf,'Stoptime',T,'SolverType','ode45');
set(conf.SolverOptions,'OutputTimes',[1:T]);
[t,pop,names]=sbiosimulate(m1);
imax=addobservable(m1,'imax','max(pop(:,2))');
sobol=sbiosobol(m1,{'beta','gamma'},'imax');
However, I am getting an error for the imax=... line which reads: "Name 'pop' in observable 'imax' does not uniquely refer to any species, parameters, compartments, or observables according to SimBiology precedence rules."
Thanks!!

Accepted Answer

Florian Augustin
Florian Augustin on 11 Nov 2021
Hi Shaun,
expressions of SimBiology observables can only reference names of components on the model, or functions on the Matlab path; referencing variables in the Matlab Workspace is not supported.
If I understand you use case correctly, you can reference "I" instead of "pop(:,2)" in the observable's expression:
addobservable(m1, "imax", "max(I)"); % <- reference I here instead of pops(:,2).
% Run Sobol analysis:
sobolResults = sbiosobol(m1, ["beta", "gamma"], "imax", "ShowWaitbar", true);
% Plot results:
bar(sobolResults);
I hope this helps.
-Florian
  6 Comments
Florian Augustin
Florian Augustin on 15 Nov 2021
Hi Shaun,
You are right, changing the (distributions of the) input parameters will generally change the Sobol indices.
Sobol indices are considered a global analysis because they reflect the sensitivity of a model's response over a whole range of input parameters (as specified by the input parameters' distributions). This is in contrast to local sensitivity analyses that only compute the sensitivity of the model response at a single point in in the parameter space.
A more specific way to describe Sobol indices would probably be "variance-based sensitivity analysis".
Sobol indices are computed with respect to the distributions of the input parameters as specified in the analysis. First order Sobol indices can be thought of as a pie chart of the total variance: the total pie is the model response's variance, and each slice (first order index) is the fraction of the variance that can be attributed to variations of a single input parameter alone. If there are interactions between parameters, then the first order Sobol indices won't explain the total variance of the model's response and we get one extra slice: "variance that is unexplained by first order Sobol indices". This is because Sobol indices are not a one-at-a-time sensitivity measure, meaning they account for interactions between input parameters.
When the distribution of the input parameters is changed, the relative effect of the sensitivity inputs changes, and as a result the Sobol indices will be different. Similarly, when input parameters are swapped out / added / removed from the analysis, then the relative effect of variations of the sensitivity inputs (and thuse the Sobol indices) changes as well. Additionally, the analysis now computes interactions between different input parameters as in the original analysis, which also changes the Sobol indices.
Best,
-Florian
Shaun Tan
Shaun Tan on 19 Nov 2021
Thanks for the explanations Florian, you've been a great help!!

Sign in to comment.

More Answers (0)

Communities

More Answers in the  SimBiology Community

Categories

Find more on Perform Sensitivity Analysis in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!