EcoSimulation — 06. Подбор параметров и поиск решения

 

Д.А. Шабанов
Конспект курса
"
Имитационное моделирование сложных биосистем
(с использованием LibreOffice Calc)
"

Имитация ограничения доступного количества ресурсов Подбор параметров и поиск решения Модели на основе клеточных автоматов
Имитационное моделирование биосистем-05 Имитационное моделирование биосистем-06 Имитационное моделирование биосистем-07

 

6. Подбор параметров и поиск решения

6.1. Подбор параметра

Мощным инструментом электронных таблиц являются средства для подбора параметров. Поясним их логику вначале на простейшем примере, а потом уже покажем их применение на примерах, где лучше проявляется их потенциал. К сожалению, с точки зрения подбора параметров Calc проигрывает Excel; для решения сложных задач нам придется перейти от одной программы к другой.

Инструмент подбора параметра Calc («Засоби/Підбір параметра...») позволяет решить следующую задачу: определить, при каком значении одной ячейки (изменяемой) значение другой ячейки (целевой) будет равно определенной величине. Естественно, это возможно в том случае, если в целевой ячейке находится формула, зависящая от значения изменяемой ячейки.

Рассмотрим такой пример. Мы имеем ряд чисел, относительно которых мы знаем (или приняли решение), какой зависимостью их следует описывать, однако не знаем значение параметра, входящего в эту зависимость. Мы хотим провести апроксимацию (поиск приближения, подходящего описания) этой зависимости. Используем, например, зависимость, описывающую экспоненциальный рост, входящую в самую первую модель, сделанную в этом курсе (пункты 2.1 – 2.4). С помощью нашей модели создадим ряд чисел, описывающих экспоненциальный рост, а потом внесем в полученный результат случайные искажения (каждую величину умножим на случайную величину в пределах 10% от самой этой величины) На рис. 6.1 вы можете увидеть, что это случайная поправка может, например, приводить к тому, что следующий ряд этой последовательности меньше предыдущего (обратите внимание на 18-й и 19-й члены ряда).

Рис. 6.1. В столбце C – значения из столбца B, измененные умножением на случайную величину в пределах ±10% (обратите внимание на строку для формул)

 

Предположим, что по значениям в столбце C нам надо установить параметры, по которым построена последовательность в столбце B (а сами значения в столбце C можно рассматривать, как измерения значений в столбце B с определенной погрешностью). Нам надо построить еще одну последовательность, которую мы будем строить в соответствии с выбранной нами закономерностью и «подгонять» под значения в столбце C. Мы сделаем это в столбце E. Заполним столбец E такой же последовательностью, как и столбец B, только вместо неизвестного нам (по столбцу C) значения r введем (в ячейку E3) значение `r – приближение к r. Мы можем задать `r каким угодно – например, равным 1, как это сделано на рис. 6.2.

Теперь мы можем сравнить имеющиеся у нас значения в столбце с C приближениями в столбце E. Для этого сравнения в столбце F высчитаем квадрат разницы между соответствующими значениями из C и E. В ячейку E4 введем сумму всех квадратов разницы между значениями из двух рядов.

Рис. 6.2. Обратите внимание на строку формул: в столбце F вычисляется квадрат разницы между соответствующими значениями из столбцов C и E

 

Теперь используем функцию «Засоби/Підбір параметра...». Появится диалог, в котором мы должны будем указать целевую ячейку с функцией, зависящей от изменяемой ячейки. Укажем, что мы хотим получить в целевой ячейке значение 0 (соответствующее идентичности двух столбцов). Изменяемая ячейка – ячейка со значением `r (рис. 6.3).

Естественно, значение в целевой ячейке не может достичь 0. Calc находит наилучшее приближение и сообщает об этом (рис. 6.4).

Рис. 6.3. Параметры диалога «Підбирання параметрів»

 

Рис. 6.4. Точное соответствие сравниваемых столбцов невозможно, но Calc подобрал неплохое приближение

 

После того, как найденное приближение будет принято, программа подставит его в ячейку для `r. Надо отметить, что приближение оказалось совсем неплохим (рис. 6.5).

Рис. 6.5. Результат подбора параметра

 

6.2. Поиск решения с несколькими параметрами: переход в Excel

Чтобы рассмотреть работу средства для поиска решения, используем данные о возрасте и длине тела 100 самок озерных лягушек, Pelophylax ridibundus, взятые из диссертационной работы Е. Е. Усовой (2016).

Таблица 6.1. Возраст (в годах) и длина тела (в мм) 100 самок озерных лягушек, Pelophylax ridibundus, (Усова, 2016)

Возраст

Длина

Возраст

Длина

Возраст

Длина

Возраст

Длина

Возраст

Длина

4

245

3

347

5

462

7

644

6

859

3

263

5

349

3

464

5

650

7

861

3

291

3

358

4

486

7

672

5

865

3

304

4

358

5

502

4

685

9

867

3

305

3

359

4

522

5

689

7

868

4

306

4

373

7

524

6

702

6

871

4

310

5

376

5

549

5

715

6

878

3

312

3

378

6

549

9

722

7

883

3

313

4

379

3

557

5

729

7

895

3

319

3

383

6

560

6

743

7

901

5

320

5

384

6

566

7

755

7

903

3

322

5

395

7

568

7

759

10

931

3

324

3

396

6

591

8

785

8

934

3

326

4

408

4

595

8

790

10

959

4

326

3

410

6

596

5

792

8

979

4

335

4

420

7

607

5

801

7

990

4

339

5

424

7

619

7

805

6

1026

4

342

5

425

6

620

8

822

7

1042

5

344

6

434

6

634

8

828

7

1055

5

345

5

447

6

642

6

852

8

1211

 

 

СТРАНИЦА В РАЗРАБОТКЕ!

 

Теперь, после того, как мы рассмотрели модели А.В. Коросова, можно применить полученный опыт к моделированию близких по структуре данных, полученных при изучении серых жаб (Bufo bufo) Иськова пруда (с. Гайдары Змиевского района Харьковской области, окрестности биостанции ХНУ).

Начиная с 2000 года каждый год батрахологи Харьковского университета собирали выборки из нерестящихся серых жаб. Каждой жабе в выборке наносили одну и ту же групповую метку, каждый год новую (мечение не проводилось в 2003 году). Начиная с 2001 года в выборках регистрировали количество жаб с метками прошлых лет. Ниже приведен фрагмент собранных данных. Он касается только самцов и охватывает промежуток времени с 2000 по 2010 гг.


  Год: 2000 2001 2002 2004 2005 2006 2007 2008 2009
  Помечено: 1257 1227 531 402 516 423 439 1055 793
Год: Выборка:                  
2001 1601 129                
2002 707 51 74              
2003 43 6 2 2            
2004 589 42 40 35            
2005 635 26 34 22 36          
2006 620 25 41 22 34 38        
2007 590 10 30 13 19 25 32      
2008 1161 6 4 3 6 20 24 14    
2009 1034 3 7 4 3 4 12 21 82  
2010 473 3 1 1 2 3 4 6 38 32

В этой таблице красным выделено количество самцов серых жаб, помеченных в каждом году. Синие цифры - выборки, исследованные на наличие в них меченых особей. Черные цифры - возвраты меток (сколько особей с метками года, подписанного красным цветом, было найдено в выборке, исследованной в году, подписанном синим цветом).

Для работы с этими данными введем следующие обозначения:

Mi – количество особей, помеченных в i-том году;

Mii+1 – количество особей, помеченных в i-том году,  в составе всей популяции в году i+1;

di+1 – доля особей в составе нерестового стада, погибших в году i+1.

Очевидно, что Mii+1 = Mi ×(1– di+1), то есть к следующему году количество особей уменьшится на долю, соответствующую смертности следующего года. К черезследующему году произойдет еще одно уменьшение численности: Mii+2 = Mi ×(1– di+1) ×(1– di+2).

Продолжим вводить обозначения:

Ni+1 – общая численность популяции в году i+1;

ni+1 – численность выборки, изученной в году i+1;

mii+1 – количество особей с метками i-того году, пойманных в году i+1 (в составе выборки численностью ni+1).

Последняя величина (mii+1) – это эмпирически определенное количество меток в выборке. Однако можно вычислить, какого количества меток следовало бы ожидать, исходя из принятых предположений. Обозначим такое математическое ожидание как m'ii+1 (обратите внимание на штрих после символа m). Как вычислить m'ii+1?

Это можно сделать на основании предположения, что доля меченых особей в выборке должна соответствовать таковой в популяции, т.е.

 m`ii+1/ni+1= Mii+1/Ni+1.

На основании этого соотношения можно установить, что 

m`ii+1 = ni+1 × Mii+1 / Ni+1,

или, с учетом ранее полученной формулы,

m`ii+1 = ni+1 × Mi ×(1 – di+1) / Ni+1.

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

N`i+1 = ni+1 × Mi ×(1 – di+1) / mii+1.

Предлагаю студентам построить модель, позволяющую определить численность самцов в данной популяции жаб и их смертность, основываясь на предположении об их постоянстве. Структура такой модели и значения вычисляемых величин показаны в файле bufo_bufo_to_the_search_for_solutions.pdf, который приложен к этой странице (внизу). Его можно не только закачать, но и открыть в окне интернет-браузера.

Как можно увидеть по pdf-файлу, в модели есть следующие блоки.

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

При построении модели желательно задать имена для тех ячеек, ссылки на которые нужно будет повторять во многих формулах. К примеру, перед ячейкой, в которой указана смертность 2001 года, укажем, что за значение в ней будет находится: напишем d2001=. Щелкнем правой кнопкой мыши на ячейке правее этой надписи, выберем "Присвоить имя", и Excel сам предложит имя d2001_ . Поскольку в модели будет использована надстройка "Поиск решения", надо, чтобы все значения переменных в модели задавались в тех ячейках, которые будет изменять надстройка. Поэтому назовем ячейку, где указывается постоянная численность, N, а ячейку для постоянной смертности — d. Затем во все ячейки, задающие численность для конкретных годов, введем формулу =N, а во все ячейки, задающие смертность для конкретных годов — формулу =d.  

После ячеек, в которых указаны численности самцов и их смертности отдельно для каждого года расположен блок с основными эмпирическими результатами. В нем указано количество помеченных особей, численности выборок, изученных в каждом году, а также количество меток каждого года в этих выборках. Этот блок организован так же, как расположенная выше таблица (за тем исключением, что в нем прописаны обозначения).

Затем следует блок, в котором расчитаны теоретически ожидаемые значения меток. К примеру, в ячейке, где находится значение m`20002001 (то есть расчетное количество особей с метками 2000 года, которых следовало бы обнаружить в выборке 2001 года) введена формула =n_2001*M2000_*(1-d2001_)/N2001_ . В ячейке для m`20022009 находится формула =n_2009*M2002_*(1-d2003_)*(1-d2004_)*(1-d2005_)*(1-d2006_)*(1-d2007_)*(1-d2008_)*(1-d2009_)/N2009_ .

В следующем блоке моделей вычислена разница между эмпирическими и теоретическими количествами вернувшихся меток. Этот блок позволяет увидеть, в каких случаях меток вернулось больше, чем можно было ожидать (им соответствуют положительные значения разностей), а в каких  — меньше (там разности отрицательные).

В следующем блока вычисляются квадраты разностей, а затем вычисляется их сумма — общая мера (Ф) несоответствия расчетных данных наблюдаемой картине.

Осталось вызвать надстройку "Поиск решения" и минимизировать значение Ф, изменяя N и d.