主要内容

sbiosobol

通过计算一阶和总阶Sobol指数执行全局敏感性分析(需要统计和机器学习工具箱)

描述

实例

sobolResults=sbiosobol(modelObj,参数个数,可观测)进行全局敏感性分析[1]在SimBiology模型上modelObj通过分解可观测关于灵敏度输入参数个数

sobolResults=sbiosobol(modelObj,情节,可观测)使用来自情节,一个SimBiology。场景对象,以执行分析。

实例

sobolResults=sbiosobol(modelObj,参数个数,可观测,名称,值)使用由一个或多个名称-值对参数指定的附加选项。

例子

全部折叠

加载肿瘤的生长模型

sbioloadproject肿瘤生长项目

获取一个变量,该变量包含要应用于模型的估计参数和剂量。

v=getvariant(m1);d=getdose(m1,“interval_dose”);

获取活动配置集并将肿瘤权重设置为响应。

c = getconfigset (m1);cs.RuntimeOptions.StatesToLog =“tumor_weight”;

模拟模型并绘制肿瘤生长曲线。

sbioplot (sbiosimulate (m1, cs, v, d));

图中包含一个轴对象。标题为“状态与时间”的axes对象包含一个类型为line的对象。该对象表示tumor_weight。

对模型进行全局敏感性分析(GSA),找出肿瘤生长敏感的模型参数。

首先,检索与肿瘤生长的药效学有关的模型参数。将模型反应定义为肿瘤重量。

modelParamNames = {“10”,“L1”,“w0”,“k1”,“k2”};输出名称=“tumor_weight”;

然后通过使用计算一阶和全阶Sobol指数来执行GSAsbiosobol.集“ShowWaitBar”真正的来显示模拟过程。默认情况下,该函数使用1000个参数样本来计算Sobol索引[1]。

rng (“默认”); sobolResults=sbiosBol(m1,modelParamName,outputName,“变异”v,“剂量”d“ShowWaitBar”,真正的)
sobolResults =具有属性的Sobol: Time: [444x1 double] SobolIndices: [5x1 struct] Variance: [444x1 table] parameters samples: [1000x5 table] Observables: {'[Tumor Growth Model]。tumor_weight `} SimulationInfo: [1x1 struct]

属性可以更改样本的数量“数字示例”名称-值对参数。该函数总共需要(输入参数个数+ 2)* NumberSamples模型模拟。

显示平均模型响应,仿真结果,以及覆盖90%仿真结果的阴影区域。

plotData (sobolResults);

图中包含一个轴对象。axis对象包含12个类型为line, patch的对象。这些对象代表模型模拟值,90.0%区域,平均值。

您可以通过指定将分位数区域调整为不同的百分比“阿尔法”对于所有模型响应的下分位数和上分位数。例如,alpha值为0.1时,在100 *α100 * (1 - α)所有模拟模型响应的分位数。

plotData (sobolResults“阿尔法”,0.1);

图中包含一个轴对象。axis对象包含12个类型为line, patch的对象。这些对象代表模型模拟,80.0%区域,平均值。

绘制一阶和全阶Sobol指数的时间历程。

h =情节(sobolResults);%调整数字大小。h、 位置(:)=[1001280800];

图中包含12个轴对象。轴对象1与标题敏感性输出[肿瘤生长模型]。Tumor_weight包含3个line类型的对象。轴对象2与标题敏感性输出[肿瘤生长模型]。Tumor_weight包含3个line类型的对象。坐标轴对象3包含3个类型为line的对象。axis对象4包含3个类型为line的对象。axis对象5包含3个类型为line的对象。轴对象6包含3个类型为line的对象。坐标轴对象7包含3个类型为line的对象。坐标轴对象8包含3个类型为line的对象。坐标轴对象9包含3个类型为line的对象。 Axes object 10 contains 3 objects of type line. Axes object 11 contains 3 objects of type line. Axes object 12 contains an object of type line.

输入参数的一阶Sobol指数给出了可单独归因于输入参数变化的总体响应方差的分数。总阶指数给出了可归因于包括输入参数变化的任何联合参数变化的总体响应方差的分数。

从Sobol指数图,参数L1w0似乎是在t = 7时应用剂量前对肿瘤重量最敏感的参数。但在施完剂量后,k1k2成为更敏感的参数,对肿瘤重量的给药后阶段贡献最大。总方差图还显示t>35时给药后阶段的方差大于给药前阶段的方差,表明k1k2可能是需要进一步研究的更重要的参数。未解释方差的分数在t = 33时显示出一定的方差,但总方差图在t = 33时显示出很小的方差,这意味着未解释方差可能不显著。未解释方差的分数计算为1 -(所有一阶Sobol指数的和),总方差计算使用风险值(响应),在那里响应为模型在每个时间点的响应。

您还可以在条形图中显示灵敏度的大小。较深的颜色意味着这些值在整个过程中出现得更频繁。

bar(sobolResults);

图中包含一个轴对象。具有标题敏感性输出的轴对象[肿瘤生长模型]。Tumor_weight包含22个patch, line类型的对象。这些对象表示一阶,总阶。

可以指定更多样本以提高Sobol索引的准确性,但模拟可能需要更长的时间才能完成。使用addsamples添加更多的样本。例如,如果指定1500个样本,函数将执行1500 *(2 +输入参数数量)模拟。

gsaMoreSamples = addsamples (gsaResults, 1500)

这个SimulationInfo属性包含计算Sobol索引的各种信息。例如,使用一组参数示例的每个仿真的模型仿真数据(SimData)存储在SimData财产的领域。这个字段是一个数组SimData对象。

sobolResults.SimulationInfo.SimData
SimBiology SimData Array: 1000-by-7 Index: Name: ModelName: dataccount: 1 - Tumor Growth Model 1 2 - Tumor Growth Model 1 3 - Tumor Growth Model 1…7000 -肿瘤生长模型1

通过检查,可以发现在计算过程中是否有模型模拟失败ValidSample领域的SimulationInfo. 在本例中,该字段显示没有失败的模拟运行。

所有(sobolResults.SimulationInfo.ValidSample)
ans =1x7逻辑阵列1 1 1 1 1

SimulationInfo。ValidSample是逻辑值的表。它的大小和SimulationInfo.SimDatA..如果ValidSample指示任何模拟失败,您可以通过从相应的列中提取信息来获取有关这些模拟运行以及用于这些运行的样本的更多信息SimulationInfo.SimDatA.假设第四列包含一个或多个失败的模拟运行。使用获取用于该模拟的模拟数据和采样值getSimulationResults

[samplesUsed, sd, validruns] = getSimulationResults (sobolResults 4);

您可以将自定义表达式添加为可观察对象,并为添加的可观察对象计算Sobol索引。例如,可以通过定义如下所示的自定义表达式来计算最大肿瘤权重的Sobol指数。

%抑制在模拟过程中发出的信息警告。警告设置=警告(“关闭”,“SimBiology: sbservices: SB_DIMANALYSISNOTDONE_MATLABFCN_UCON”);%添加可观察表达式。sobolObs = addobservable (sobolResults,“最大tumor_weight”,“马克斯(tumor_weight)”,“单位”,“克”);

绘制显示90%分位数区域的计算仿真结果。

h2 = plotData (sobolObs);h1 . position (:) = [100 100 1280 800];

图中包含2个轴对象。坐标轴对象1包含12个类型为line, patch的对象。这些对象代表模型模拟值,90.0%区域,平均值。axis对象2包含12个类型为line, patch的对象。这些对象代表模型模拟值,90.0%区域,平均值。

您还可以通过指定可观察对象的名称来删除它。

gsaNoObs = removeobservable (sobolObs,“最大tumor_weight”);

恢复警告设置。

警告(warnSettings);

输入参数

全部折叠

SimBiology模型,指定为SimBiology模型对象

模型参数、物种或单元的名称,指定为字符向量、字符串、字符串向量或字符向量的单元数组。

例子:[“k1”、“k2”]

数据类型:字符|字符串|细胞

模型响应,指定为字符向量、字符串、字符串向量或字符向量的单元数组。指定物种、参数、隔间或的名称可观测

例子:[" tumor_growth "]

数据类型:字符|字符串|细胞

图纸样品的来源,指定为SimBiology。场景对象。

  • 必须按元素组合对象的项。

  • 条目必须是独立的随机变量。如果有多个条目,它们必须是不相关的。

  • SamplingMethod任何条目都不能是连系动词

名称值参数

指定可选的逗号分隔的字符对名称,值参数。的名字是参数名和价值为对应值。的名字必须出现在引号内。可以以任意顺序指定多个名称和值对参数Name1, Value1,…,的家

例子:sobolResults=sbiosBol(modelObj,参数,可观测值,'ShowWaitbar',true)指定显示模拟进度条。

参数边界,指定为具有两列的数字矩阵。第一列包含下限,第二列包含上限。行数必须等于中的参数数参数个数

如果参数具有非零值,则默认边界为该值的±10%。如果参数值为零,则默认边界为[0 1]

例子:(0.5 - 5)

数据类型:

在模型模拟期间使用的剂量,指定为计划剂量RepeatDose物体或剂量物体的载体。

在模型模拟之前应用的变量,指定为变量对象或变量对象的向量。

当为属性值指定多个具有重复规范的变体时,将在模拟期间使用该属性值在变体数组中的最后一次出现。

要计算Sobol索引的样本数,指定为逗号分隔对,由“数字示例”一个正整数。功能要求(数量的输入参数个数+ 2) *数字样本模型模拟,计算一阶和全阶Sobol指数。

数据类型:

用于绘制样本的概率分布,指定为概率。ProbabilityDistribution这些对象的对象或向量。指定标量概率。ProbabilityDistribution或者是长度向量N,在那里N是输入参数的个数。可以使用以下方法创建分布对象,以从各种分布(如均匀分布、正态分布或对数正态分布)中进行抽样makedist(统计和机器学习工具箱)

如果指定标量概率。ProbabilityDistribution对象,并且有多个输入参数,sbiosobol使用相同的分布对象为每个参数绘制采样。

你不能同时指定这个参数界限

SimBiology。场景对象是输入。

方法生成参数样本,该参数样本被指定为字符向量或字符串。有效的选项是:

  • “Sobol”—使用低差异Sobol序列生成样本。

  • “荷”-使用低差异Halton序列生成样本。

  • “韩”-使用低差异的拉丁超立方体样本。

  • “随机”—使用均匀分布的随机样本。

数据类型:字符|字符串

采样方法的选项,指定为标量结构体。根据抽样方法的不同,选择也不同:sobol,哈尔顿,或

对于sobol哈尔顿的每个name-value参数指定结构的每个字段名和值索波尔塞特(统计和机器学习工具箱)haltonset(统计和机器学习工具箱)函数。SimBiology使用默认值1.跳过两个方法的参数。对于所有其他名称-值参数,软件使用相同的默认值索波尔塞特haltonset。例如,为飞跃跳过具有以下非默认值的选项。

s1.Leap=50;s1.Skip=0;

对于,有三个采样器支持不同的采样选项。金宝app

  • 如果你指定一个协方差矩阵,SimBiology使用lhsnorm(统计和机器学习工具箱)抽样。SamplingOptions参数不允许。

  • 否则,使用字段名使用设计选择取样器。

    • 如果值为真正的, SimBiology使用lhsdesign(统计和机器学习工具箱).的名称-值参数lhsdesign指定字段名称和值。

    • 如果值为错误的(默认),SimBiology使用一个不可配置的拉丁超立方体采样器lhsdesign.该采样器不需要统计学和机器学习工具箱™。SamplingOptions不能包含任何其他选项,除非使用设计

例如,设置一个要使用的结构lhsdesign标准迭代选项。

s2。UseLhsdesign = true;s2。标准=“相关”; s2.迭代次数=10次;

SimBiology。场景对象是输入。

模拟停止时间,指定为非负标量。如果未指定结束时刻也不OutputTimes,该函数使用模型的活动配置集中的停止时间。不能同时指定这两个时间结束时刻OutputTimes

数据类型:

仿真输出次数,指定为逗号分隔对组成“OutputTimes”和数字向量。该函数计算这些输出时间点的Sobol索引。不能同时指定两者结束时刻OutputTimes.默认情况下,该函数使用第一个模型模拟的输出次数。

例子:[0 1 2 3.5 4 5 5.5]

数据类型:

标志以并行运行模型模拟,指定为真正的错误的.当值为真正的和Parallel Computing Toolbox™可用时,该函数并行运行模拟。

数据类型:必然的

打开模型加速的标志,指定为真正的错误的

数据类型:必然的

将模型响应插值到一组公共输出时间的方法,指定为字符向量或字符串。下面是有效的选项。

  • “interp1q”——使用interp1q函数。

  • 使用interp1函数,指定以下方法之一:

    • “最近的”

    • “线性”

    • “样条曲线”

    • “pchip”

    • “v5cubic”

  • “zoh”-指定零阶保持。

数据类型:字符|字符串

通过显示进度条来显示模型模拟的进度的标志,该进度条指定为逗号分隔对,由“ShowWaitbar”真正的错误的.缺省情况下,不显示等待条。

数据类型:必然的

输出参数

全部折叠

结果包含第一阶和总阶Sobol索引,返回为SimBiology.gsa.Sobol对象。该对象还包含用于计算Sobol指数的参数样本值和模型模拟数据。

结果对象可以包含大量的模拟数据(SimData)。对象的大小超过(1 +观测值个数)*输出时间点个数*(2 +参数个数)*样本个数* 8例如,如果您有一个可观察的、500个输出时间点、8个参数和100000个样本,则对象大小为(1 + 1) * 500 * (2 + 8) * 100000 * 8bytes = 8gb。如果你需要保存这样大的对象,请使用下面的语法:

保存(variableName文件名,“-v7.3”);
有关详细信息,请参阅MAT文件版本

更多关于

全部折叠

计算Sobol指数的Saltelli方法

sbiosobol实现Saltelli方法[1]计算Sobol指数。

考虑一个SIM生物学模型响应Y用数学模型表示的 Y = F ( X 1. , X 2. , X 3. , ... , X K ) ,在那里X是一个模型参数和我= 1,…,k

一阶Sobol指数(s)给出总体响应方差的分数V (Y)这可以归因于X一个人。s定义如下。

s = v X ( E X ( Y | X ) ) v ( Y )

Sobol总顺序索引(s“透明国际”)给出总体响应方差的分数V (Y)这可以归因于任何关节参数的变化包括的变化Xs“透明国际”定义如下。

s T = 1. v X ( E X ( Y | X ) ) v ( Y ) = E X ( v X ( Y | X ) ) v ( Y )

计算…的个别值Y对应样品的参数X1.,X1.、……XK考虑两个独立的采样矩阵A.B

A. = ( X 11 X 12 ... X 1. K X 21 X 22 ... X 2. K ... ... ... ... X N 1. X N 2. ... X N K )

B = ( X 11 ' X 12 ' ... X 1. K ' X 21 ' X 22 ' ... X 2. K ' ... ... ... ... X N 1. ' X N 2. ' ... X N K ' )

N为样本量。矩阵的每一行A.B对应一个参数样本集,是模型参数值的单一实现。

预估ss“透明国际”是从模型模拟结果中使用矩阵的样本值得到的A.,B, A. B ,它是一个所有列都来自的矩阵A.除了第th列,它来自B对于I = 1, 2,…,参数

A. B = ( X 11 X 12 ... X 1. ' ... X 1. K X 21 X 22 ... X 2. ' ... X 2. K ... ... ... ... ... ... X N 1. X N 2. ... X N ' ... X N K )

一阶Sobol指数和全阶Sobol指数的近似公式如下:

s ^ = 1. N J = 1. N F ( B ) J ( F ( A. B ) J F ( A. ) J ) v ( Y )

s ^ T = 1. 2. N J = 1. N ( F ( A. ) J F ( A. B ) J ) 2. v ( Y )

f(A.),f(B), F ( A. B ) J 模型模拟结果是否使用矩阵中的参数样本值A.,B, A. B

矩阵A.对应于ParameterSamples属性的Sobol结果对象(resultsObj。ParameterSamples).矩阵B对应于金宝app支持样本财产(resultsObj.SimulationInfo.金宝appSupportSamples).

这个 A. B 矩阵存储在SimData的结构SimulationInfo财产(resultsObj.SimulationInfo.SimData).大小SimulationInfo.SimDatA.NumberSamples——- - - - - -参数个数+ 2,在那里NumberSamples样品的数量是多少参数是输入参数的个数。列数是2 +参数个数因为第一列SimulationInfo.SimDatA.包含使用样本矩阵的模型仿真结果A..第二列包含使用金宝app支持样本,这是另一个样本矩阵B。其余列包含使用的模拟结果 A. B 1. , A. B 2. 、…… A. B 、…… A. B P A. R A. M s .看到getSimulationResults检索指定的模型仿真结果和样本th指数( A. B )SimulationInfo.SimDatA.数组中。

参考文献

[1]萨特利,安德里亚,保拉·安诺尼,伊万诺·阿兹尼,弗朗西斯卡·坎波隆戈,马可·拉托和斯特凡诺·塔兰托拉。“基于方差的模型输出灵敏度分析”。总敏感性指数的设计和估计。计算机物理通信181年,没有。2(2010年2月):259-70。https://doi.org/10.1016/j.cpc.2009.09.018。

介绍了R2020a