sbiosobol

Выполните глобальный анализ чувствительности путем вычисления сначала - и общий порядок индексы Sobol (требует Statistics and Machine Learning Toolbox),

Описание

пример

sobolResults = sbiosobol(modelObj,params,observables) выполняет глобальный анализ чувствительности [1] на модели SimBiology modelObj путем разложения отклонений observables относительно чувствительности вводит params.

пример

sobolResults = sbiosobol(modelObj,params,observables,Name,Value) дополнительные опции использования заданы одним или несколькими аргументами пары "имя-значение".

Примеры

свернуть все

Загрузите модель роста опухоли.

sbioloadproject tumor_growth_vpop_sa.sbproj

Заставьте вариант предполагаемыми параметрами и дозой применяться к модели.

v = getvariant(m1);
d = getdose(m1,'interval_dose');

Получите активный configset и установите вес опухоли как ответ.

cs = getconfigset(m1);
cs.RuntimeOptions.StatesToLog = 'tumor_weight';

Симулируйте модель и постройте профиль роста опухоли.

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

Figure contains an axes. The axes with title States versus Time contains an object of type line. This object represents tumor_weight.

Выполните глобальный анализ чувствительности (GSA) на модели, чтобы найти параметры модели, к которым рост опухоли чувствителен.

Во-первых, получите параметры модели интереса, которые вовлечены в фармакодинамику роста опухоли. Задайте ответ модели как вес опухоли.

modelParamNames = {'L0','L1','w0','k1','k2'};
outputName = 'tumor_weight';

Затем выполните GSA путем вычисления первого - и общий порядок индексы Sobol с помощью sbiosobol. Установите 'ShowWaitBar' к true показать прогресс симуляции. По умолчанию функция использует 1 000 выборок параметра, чтобы вычислить индексы Sobol [1].

sobolResults = sbiosobol(m1,modelParamNames,outputName,'Variants',v,'Doses',d,'ShowWaitBar',true)
sobolResults = 
  Sobol with properties:

                Time: [444x1 double]
        SobolIndices: [5x1 struct]
            Variance: [444x1 table]
         Observables: {'[Tumor Growth Model].tumor_weight'}
    ParameterSamples: [1000x5 table]
      SimulationInfo: [1x1 struct]

Можно изменить количество выборок путем определения 'NumberSamples' аргумент пары "имя-значение". Функция требует в общей сложности (number of input parameters + 2) * NumberSamples симуляции модели.

Покажите средний ответ модели, результаты симуляции и теневую область, покрывающую 90% результатов симуляции.

plotData(sobolResults);

Figure contains an axes. The axes contains 12 objects of type line, patch. These objects represent model simulation, 90.0% region, mean value.

Можно настроить область квантиля к различному проценту путем определения 'Alphas' для более низких и верхних квантилей всех ответов модели. Например, альфа-значение 0,1 графиков теневая область между 100 * alpha и 100 * (1 - alpha) квантили всех симулированных ответов модели.

plotData(sobolResults,'Alphas',0.1);

Figure contains an axes. The axes contains 12 objects of type line, patch. These objects represent model simulation, 80.0% region, mean value.

Постройте курс времени первого - и общий порядок индексы Sobol.

h = plot(sobolResults);
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Figure contains 12 axes. Axes 1 with title [Tumor Growth Model].tumor_weight contains 3 objects of type line. Axes 2 with title [Tumor Growth Model].tumor_weight contains 3 objects of type line. Axes 3 contains 3 objects of type line. Axes 4 contains 3 objects of type line. Axes 5 contains 3 objects of type line. Axes 6 contains 3 objects of type line. Axes 7 contains 3 objects of type line. Axes 8 contains 3 objects of type line. Axes 9 contains 3 objects of type line. Axes 10 contains 3 objects of type line. Axes 11 contains 3 objects of type line. Axes 12 contains an object of type line.

Индекс Sobol первого порядка входного параметра дает часть полного отклонения ответа, которое может быть приписано изменениям одного только входного параметра. Индекс общего порядка дает часть полного отклонения ответа, которое может быть приписано любым объединенным изменениям параметра, которые включают изменения входного параметра.

Из графиков индексов Sobol, параметры L1 и w0 кажется, самые чувствительные параметры к весу опухоли, прежде чем доза была применена в t = 7. Но после того, как доза применяется, k1 и k2 станьте более чувствительными параметрами и способствуйте больше всего дозируемому после того, как этапу веса опухоли. Общий график отклонения также показывает большее отклонение для этапа после-того,-как-дозы в t> 35, чем для этапа прежде-чем-дозы роста опухоли, указывая на тот k1 и k2 могут быть более важные параметры, чтобы заняться расследованиями далее. Часть необъясненного отклонения показывает некоторое отклонение приблизительно в t = 33, но общий график отклонения показывает мало отклонения в t = 33, означая, что необъясненное отклонение могло быть незначительным. Часть необъясненного отклонения вычисляется как 1 - (сумма всех индексов Sobol первого порядка), и общее отклонение вычисляется с помощью var(response), где response ответ модели в каждом моменте времени.

Можно также отобразить величины чувствительности в столбиковой диаграмме.

bar(sobolResults)

Figure contains an axes. The axes with title [Tumor Growth Model].tumor_weight contains 22 objects of type patch, line. These objects represent first order, total order.

Можно задать больше выборок, чтобы увеличить точность индексов Sobol, но симуляция может занять больше времени, чтобы закончиться. Используйте addsamples добавить больше выборок. Например, если вы задаете 1 500 выборок, функция выполняет 1500 * (2 + number of input parameters) симуляции.

gsaMoreSamples = addsamples(gsaResults,1500)

Свойство SimulationInfo объекта результата содержит различную информацию для вычисления индексов Sobol. Например, данные о симуляции модели (SimData) для каждой симуляции с помощью набора выборок параметра хранятся в SimData поле свойства. Это поле является массивом SimData объекты.

sobolResults.SimulationInfo.SimData
 
   SimBiology SimData Array : 1000-by-7
 
   Index:    Name:         ModelName:         DataCount: 
   1           -           Tumor Growth Model 1          
   2           -           Tumor Growth Model 1          
   3           -           Tumor Growth Model 1          
   ...                                                   
   7000        -           Tumor Growth Model 1          
 

Можно узнать, перестала ли какая-либо симуляция модели работать во время расчета путем проверки ValidSample поле SimulationInfo. В этом примере поле не показывает не пройдено запусков симуляции.

all(sobolResults.SimulationInfo.ValidSample)
ans = 1x7 logical array

   1   1   1   1   1   1   1

SimulationInfo.ValidSample таблица логических значений. Это имеет тот же размер как SimulationInfo.SimData. Если ValidSample указывает, что любые симуляции перестали работать, можно получить больше информации о тех запусках симуляции и выборках, используемых для тех запусков путем извлечения информации из соответствующего столбца SimulationInfo.SimData. Предположим, что четвертый столбец содержит один или несколько не пройдено запусков симуляции. Получите данные моделирования и демонстрационные значения, используемые для той симуляции с помощью getSimulationResults.

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

Можно добавить пользовательские выражения как observables и вычислить индексы Sobol для добавленного observables. Например, можно вычислить индексы Sobol для максимального веса опухоли путем определения пользовательского выражения можно следующим образом.

% Suppress an information warning that is issued during simulation.
warnSettings = warning('off', 'SimBiology:sbservices:SB_DIMANALYSISNOTDONE_MATLABFCN_UCON');
% Add the observable expression.
sobolObs = addobservable(sobolResults,'Maximum tumor_weight','max(tumor_weight)','Units','gram');

Постройте вычисленные результаты симуляции, показывающие 90%-ю область квантиля.

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

Figure contains 2 axes. Axes 1 contains 12 objects of type line, patch. These objects represent model simulation, 90.0% region, mean value. Axes 2 contains 12 objects of type line, patch. These objects represent model simulation, 90.0% region, mean value.

Можно также удалить заметное путем определения его имени.

gsaNoObs = removeobservable(sobolObs,'Maximum tumor_weight');

Восстановите настройки предупреждения.

warning(warnSettings);

Входные параметры

свернуть все

Модель SimBiology в виде SimBiology model object.

Имена параметров модели, разновидностей или отсеков в виде вектора символов, строки, представляют в виде строки вектор или массив ячеек из символьных векторов.

Пример: ["k1","k2"]

Типы данных: char | string | cell

Ответы модели в виде вектора символов, строки, представляют в виде строки вектор или массив ячеек из символьных векторов. Задайте имена разновидностей, параметров, отсеков или observables.

Пример: ["tumor_growth"]

Типы данных: char | string | cell

Аргументы в виде пар имя-значение

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: sobolResults = sbiosobol(modelObj,params,observables,'ShowWaitbar',true) задает, чтобы показать индикатор выполнения симуляции.

Параметр ограничивает в виде разделенной запятой пары, состоящей из 'Bounds' и числовая матрица с двумя столбцами. Первый столбец содержит нижние границы, и второй столбец содержит верхние границы. Количество строк должно быть равно количеству параметров в params.

Если параметр имеет ненулевое значение, границы по умолчанию составляют ±10% значения. Если значение параметров является нулем, границами по умолчанию является [0 1].

Пример: [0.5 5]

Типы данных: double

Дозы, чтобы использовать во время симуляций модели в виде разделенной запятой пары, состоящей из 'Doses' и ScheduleDose или RepeatDose возразите или вектор из объектов дозы.

Варианты, чтобы применяться перед симуляциями модели в виде разделенной запятой пары, состоящей из 'Variants' и различный объект или вектор из различных объектов.

Когда существует несколько вариантов и если существуют дублирующиеся технические требования для значения свойства, последнее вхождение для значения свойства в массиве вариантов используется в процессе моделирования.

Количество выборок, чтобы вычислить индексы Sobol в виде разделенной запятой пары, состоящей из 'NumberSamples' и положительное целое число. Функция требует (количество входа params + 2) * NumberSamples симуляции модели, чтобы вычислить первое - и общий порядок индексы Sobol.

Типы данных: double

Метод, чтобы сгенерировать выборки параметра в виде разделенной запятой пары, состоящей из 'SamplingMethod' и вектор символов или строка. Допустимые опции:

  • 'Sobol' — Используйте низкое несоответствие последовательность Sobol, чтобы сгенерировать выборки.

  • 'Halton' — Используйте низкое несоответствие последовательность Холтона, чтобы сгенерировать выборки.

  • 'lhs' — Используйте латинские выборки гиперкуба низкого несоответствия.

  • 'random' — Используйте равномерно распределенные случайные выборки.

Типы данных: char | string

Время остановки симуляции в виде разделенной запятой пары, состоящей из 'StopTime' и неотрицательный скаляр. Если вы не задаете никакой StopTime ни OutputTimes, функция использует время остановки от активной конфигурации модели модели. Вы не можете задать оба StopTime и OutputTimes.

Типы данных: double

Simulation времена выхода в виде разделенной запятой пары, состоящей из 'OutputTimes' и числовой вектор. Функция вычисляет индексы Sobol в этих выходных моментах времени. Вы не можете задать оба StopTime и OutputTimes. По умолчанию функция использует выходные времена первой симуляции модели.

Пример: [0 1 2 3.5 4 5 5.5]

Типы данных: double

Отметьте, чтобы запустить симуляции модели параллельно в виде разделенной запятой пары, состоящей из 'UseParallel' и true или false.

Типы данных: логический

Отметьте, чтобы включить ускорение модели в виде разделенной запятой пары, состоящей из 'Accelerate' и true или false.

Типы данных: логический

Метод для интерполяции симуляций модели к новым выходным временам в виде разделенной запятой пары, состоящей из 'InterpolationMethod' и вектор символов или строка. Допустимые опции следуют.

  • "interp1q" — Используйте interp1q функция.

  • Используйте interp1 функция путем определения одного из следующих методов:

    • "nearest"

    • "linear"

    • "spline"

    • "pchip"

    • "v5cubic"

  • "zoh" — Укажите, что нулевой порядок содержит.

Типы данных: char | string

Отметьте, чтобы показать прогресс симуляций модели путем отображения индикатора выполнения в виде разделенной запятой пары, состоящей из 'ShowWaitbar' и true или false. По умолчанию никакие не ожидают, панель отображена.

Типы данных: логический

Выходные аргументы

свернуть все

Результаты, содержащие первое - и общий порядок индексы Sobol, возвращенные как SimBiology.gsa.Sobol объект. Объект также содержит демонстрационные значения параметра, и данные о симуляции модели использовались для расчета индексов Sobol.

Объект результатов может содержать существенное количество данных моделирования (SimData). Размер объекта превышает (1 + number of observables) * number of output time points * (2 + number of parameters) * number of samples * 8 байты. Например, если у вас есть одно заметное, 500 выходных моментов времени, 8 параметров и 100 000 выборок, размером объекта является (1 + 1) * 500 * (2 + 8) * 100000 * 8 байты = 8 Гбайт. Если необходимо сохранить такие большие объекты, используйте этот синтаксис:

save(fileName,variableName,'-v7.3');
Для получения дополнительной информации смотрите версию MAT-файла.

Больше о

свернуть все

Метод Сальтелли, чтобы вычислить индексы Sobol

sbiosobol реализует метод Сальтелли [1], чтобы вычислить индексы Sobol.

Считайте ответ модели SimBiology Y описанный как математическая модель Y=f(X1,X2,X3,...,Xk), где Xi является параметром модели и i = 1,…,k.

Индекс Sobol первого порядка (Si) дает часть полного отклонения ответа V(Y) это может быть приписано изменениям одного только Xi. Si определяется следующим образом.

Si=VXi(EXi(Y|Xi))V(Y)

Общий порядок индекс Sobol (STi) дает часть полного отклонения ответа V(Y) это может быть приписано любым объединенным изменениям параметра, которые включают изменения Xi. STi определяется следующим образом.

STi=1VXi(EXi(Y|Xi))V(Y)=EXi(VXi(Y|Xi))V(Y)

Чтобы вычислить отдельные значения для Y, соответствующего выборкам параметров X1, X1, …, Xk, считает две независимых матрицы выборки A и B.

A=(X11X12...X1kX21X22...X2k............Xn1Xn2...Xnk)

B=(X11'X12'...X1k'X21'X22'...X2k'............Xn1'Xn2'...Xnk')

n является объемом выборки. Каждая строка матриц, A и B соответствуют одному демонстрационному набору параметра, который является одной реализацией значений параметра модели.

Оценки для Si и STi получены из результатов симуляции модели с помощью демонстрационных значений из матриц A, B, и ABi, который является матрицей, где все столбцы являются от A кроме i th столбцом, который является от B для i = 1, 2, …, params.

ABi=(X11X12...X1i'...X1kX21X22...X2i'...X2k..................Xn1Xn2...Xni'...Xnk)

Формулы, чтобы аппроксимировать первое - и общий порядок индексы Sobol следующие.

S^i=1nj=1nf(B)j(f(ABi)jf(A)j)V(Y)

S^Ti=12nj=1n(f(A)jf(ABi)j)2V(Y)

f(A), f(B), и f(ABi)j результаты симуляции модели с помощью демонстрационных значений параметра из матриц A, B, и ABi.

Матричный A соответствует ParameterSamples свойство Sobol заканчивается объект (resultsObj.ParameterSamples). Матричный B соответствует SupportSamples свойство (resultsObj.SimulationInfo.SupportSamples).

ABi матрицы хранятся в SimData структура SimulationInfo свойство (resultsObj.SimulationInfo.SimData). Размер SimulationInfo.SimData NumberSamples-by-params + 2, где NumberSamples является количеством выборок, и param является количеством входных параметров. Количеством столбцов является 2 + params потому что первый столбец SimulationInfo.SimData содержит результаты симуляции модели с помощью демонстрационного матричного A. Второй столбец содержит результаты симуляции с помощью SupportSamples, который является другим демонстрационным матричным B. Остальная часть столбцов содержит использование результатов симуляции AB1, AB2, …, ABi, …, ABparams. Смотрите getSimulationResults получать результаты симуляции модели и выборки для заданного i th индекс (ABi) от SimulationInfo.SimData массив.

Ссылки

[1] Сальтелли, Андреа, Паола Аннони, Ivano Azzini, Франческа Камполонго, Марко Ратто и Стефано Тарантола. “Основанный на отклонении Анализ чувствительности Выхода Модели. Проект и Средство оценки для Общего индекса Чувствительности”. Computer Physics Communications 181, № 2 (февраль 2010): 259–70. https://doi.org/10.1016/j.cpc.2009.09.018.

Введенный в R2020a