Подгонка смоделированных и экспериментальных точек данных с Python

Question

Я написал некоторый код, который выполняет моделирование по методу Монте-Карло и выдает кривые зависимости интенсивности сигнала от времени.Форма такой кривой зависит от различных параметров, два из которых мой сотрудник хочет определить по «реальной версии» эксперимента, который я моделирую.

Мы готовы сравнить ее экспериментальные данные с моими смоделированными кривыми, но теперь я застрял, поскольку я еще не смог выполнить какую-либо подгонку (пока я заменил экспериментальные данные смоделированными данными с шумом длятестирование).Я попытался использовать scipy.optimize.leastsq, который завершается с кодом 2 (в соответствии с документацией это означает, что подгонка прошла успешно), но в основном он просто возвращает значения (не совсем те же, но близкие), которые я ввел в качестве начального предположения,независимо от того, насколько близко или далеко они были от истинных ценностей.Если он сообщает о других значениях, результирующая кривая все еще значительно отличается от истинной.

Другое наблюдение состоит в том, что infodict['nfev'] неизменно содержит

The relative error between two consecutive iterates is at most 0.000000

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

Кто-нибудь знает, что я могу делать неправильно, или какую функцию подгонки я мог бы использовать вместо leastsq?Если так: большое спасибо заранее!

РЕДАКТИРОВАТЬ

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

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

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

РЕДАКТИРОВАТЬ 2

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

РЕДАКТИРОВАТЬ 3

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

Answer 1

Для подгонки ваших данных к произвольной функции вам обычно нужен алгоритм Левенберга – Марквардта , и именно это использует scipy.optimize.leastsq, поэтому вы, скорее всего, используете правильную функцию.

Посмотрите учебник в разделе 5.4 этой страницы и посмотрите, поможет ли он.

Также возможно, что ваша основная функция трудна для подгонки (что это за функция?), Но, похоже, у вас могут быть другие проблемы.

Кроме того, как бы ни был хорош StackOverflow, вы, вероятно, получите гораздо более квалифицированные ответы на скупые вопросы, отправив их непосредственно в список рассылки Scipy-User с некоторым примером кода и более подробной информацией.

Answer 2

Не совсем ответ, но есть также PyMinuit, чтобы попробовать.

http://code.google.com/p/pyminuit/

То, что вы хотите сделать, - это преобразовать ваши pdf и точки данных в chi ^ 2 или -ln (правдоподобие) или метрику по вашему выбору и использовать PyMinuit для минимизации этой метрики. Его можно настроить как очень подробный, чтобы вы могли узнать, где что-то пошло не так (если оно пошло не так).

Answer 3

Если вы не знаете ожидаемую функциональную форму глобального объекта, но может предсказать ожидаемое значение для «следующей» точки, учитывая текущее состояние системы, которое вы могли бы рассмотреть, используя Фильтр Калмана (да, «фильтр» звучит глупо в подходящем контексте, но название является историческим и не может быть легко изменено сейчас).

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

1. Определить пространство представления
2. Выразите ваши данные (смоделированные или экспериментальные) в пространстве представления.
3. Определить процедуру обновления, которая получает «следующее» представление от заданного
4. Извлечение требуемой кривой из представления серии, возвращаемого установщиком

Кажется, по крайней мере один существующий пакет Python поддерживает это (обратите внимание, что интерфейс здесь отличается от того, к которому я привык, и я не могу дать много советов по его использованию ).

Answer 4

Поскольку у вас есть только два параметра, вам нужно просто выполнить поиск по сетке.

results = {}
for p0 in parameter_space_for_first_parameter:
     for p1 in parameter_space_for_second_parameter:
           results[p0,p1] = compare(p0,p1)

Если вы можете позволить себе вычисления, compare должен выполнить несколько прогонов (с разными инициализациями) и вычислить среднее и стандартное отклонения. Вы можете попробовать использовать мой пакет jug для управления своими вычислениями (он был разработан именно для такого рода вещей).

Наконец, выведите результаты, посмотрите на минимумы (их может быть несколько). Это «глупый» метод, но он будет работать во многих ситуациях, когда другие методы застревают.

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