Formatting for 'addobservables'
显示 更早的评论
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!!
采纳的回答
Florian Augustin
2021-11-11
1 个投票
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 个评论
Shaun Tan
2021-11-12
Shaun Tan
2021-11-12
Florian Augustin
2021-11-12
Hi Shaun,
You are very much on the right track with the observable.
If I understand the use-case correctly, the observable should return the time at which the max. I is attained. You can do this by indexing into the time variable, which is a reserved token in SimBiology that gives you access to simulation times of the ODE solver.
When evaluating expressions of observables, SimBiology first tries to determine whether the result is a scalar value or a vector. I therefore recommend to write observables' expressions in a way that they either always return a scalar value or always a vector of values, for any value of I. In your case, find(I == max(I)), returns a scalar if the max. is unique, but it can also evaluate to a vector if the max. is not unique.
If you know that the max. is unique in all your model simulations, you can force the index to be unique by choosing the first occurance as follows:
addobservable(m1, "iday", "time(find(I == max(I), 1))");
If you are unsure and your model could attain the max. value at several different time points and you want to be alerted if this happens, you can write a helper function (an example, getTimeOfMaxValue.m, is attached) to guard against this case:
addobservable(m1, "iday", "getTimeOfMaxValue(time, I)");
Best,
Florian
Shaun Tan
2021-11-13
Florian Augustin
2021-11-15
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
2021-11-19
更多回答（0 个）
社区
更多回答在 SimBiology Community
类别
在 帮助中心 和 File Exchange 中查找有关 Scan Parameter Ranges 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
