аппроксимация экспоненциального распада без начального предположения

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

>>> import numpy as np
>>> from scipy.optimize import curve_fit
>>> def func(x, a, b, c):
...     return a*np.exp(-b*x) + c

>>> x = np.linspace(0,4,50)
>>> y = func(x, 2.5, 1.3, 0.5)
>>> yn = y + 0.2*np.random.normal(size=len(x))

>>> popt, pcov = curve_fit(func, x, yn)

Подгоненные параметры будут отличаться из-за добавленного случайного шума, но я получил 2.47990495, 1.40709306, 0.53753635 как a, b и c, так что это не так уж плохо с шумом там. Если я подхожу к y вместо yn, я получаю точные значения a, b и c.

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

Вот краткий обзор

http://www.statsci.org/other/prony.html

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

Оценка Прони требует, чтобы смещение было известно, но если вы идете «достаточно далеко» в своем затухании, у вас есть разумная оценка смещения, так что вы можете просто сдвинуть данные, чтобы сместить их на 0. В любом Скорость, оценка Прони - это просто способ получить разумное начальное предположение для других подходящих процедур.

Я не знаю python, но я знаю простой способ не итеративной оценки коэффициентов экспоненциального затухания со смещением, учитывая три точки данных с фиксированной разницей в их независимых координатах. Ваши точки данных имеют фиксированную разницу в своих независимых координатах (ваши значения x расположены с интервалом 60), поэтому мой метод может быть применен к ним. Вы, конечно, можете перевести математику в Python.

Предположим

y = A + B*exp(-c*x) = A + B*C^x

где C = exp(-c)

Учитывая y_0, y_1, y_2, для x = 0, 1, 2, мы решаем

y_0 = A + B
y_1 = A + B*C
y_2 = A + B*C^2

найти A, B, C следующим образом:

A = (y_0*y_2 - y_1^2)/(y_0 + y_2 - 2*y_1)
B = (y_1 - y_0)^2/(y_0 + y_2 - 2*y_1)
C = (y_2 - y_1)/(y_1 - y_0)

Соответствующая экспонента проходит точно через три точки (0, y_0), (1, y_1) и (2, y_2). Если ваши точки данных находятся не в координатах x 0, 1, 2, а в точках k, k + s и k + 2 * s, то

y = A′ + B′*C′^(k + s*x) = A′ + B′*C′^k*(C′^s)^x = A + B*C^x

, так что вы можете использовать вышеуказанные формулы, чтобы найти A , B, C, а затем вычислите

A′ = A
C′ = C^(1/s)
B′ = B/(C′^k)

Полученные коэффициенты очень чувствительны к ошибкам в координатах y, которые могут привести к большим ошибкам, если вы экстраполируете за пределы диапазона, определенного тремя используемыми точками данных, поэтому Лучше всего рассчитывать A, B, C из трех точек данных, которые находятся как можно дальше друг от друга (при этом между ними все еще остается фиксированное расстояние).

Ваш набор данных имеет 10 равноудаленных точек данных. Давайте выберем три точки данных (110, 2391), (350, 786), (590, 263) для использования - они имеют максимально возможное фиксированное расстояние (240) в независимой координате. Итак, y_0 = 2391, y_1 = 786, y_2 = 263, k = 110, s = 240. Тогда A = 10.20055, B = 2380.799, C = 0.3258567, A ′ = 10.20055, B ′ = 3980.329, C ′ = 0.9953388. Экспонента равна

y = 10.20055 + 3980.329*0.9953388^x = 10.20055 + 3980.329*exp(-0.004672073*x)

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

Формула для вычисления A такая же, как и при преобразовании Шенкса ( http://en.wikipedia.org/wiki/Shanks_transformation ).

Python-реализация решения @ JJacquelin. Мне нужно было приблизительное решение без решения без начальных догадок, поэтому ответ @ JJacquelin был действительно полезным. Первоначальный вопрос был задан как запрос питона на numpy / scipy. Я взял хороший чистый код R @ johanvdw и реорганизовал его как python / numpy. Надеюсь, кому-нибудь пригодится: https://gist.github.com/friendtogeoff/00b89fa8d9acc1b2bdf3bdb675178a29

import numpy as np

"""
compute an exponential decay fit to two vectors of x and y data
result is in form y = a + b * exp(c*x).
ref. https://gist.github.com/johanvdw/443a820a7f4ffa7e9f8997481d7ca8b3
"""
def exp_est(x,y):
    n = np.size(x)
    # sort the data into ascending x order
    y = y[np.argsort(x)]
    x = x[np.argsort(x)]

    Sk = np.zeros(n)

    for n in range(1,n):
        Sk[n] = Sk[n-1] + (y[n] + y[n-1])*(x[n]-x[n-1])/2
    dx = x - x[0]
    dy = y - y[0]

    m1 = np.matrix([[np.sum(dx**2), np.sum(dx*Sk)],
                    [np.sum(dx*Sk), np.sum(Sk**2)]])
    m2 = np.matrix([np.sum(dx*dy), np.sum(dy*Sk)])

    [d, c] = (m1.I * m2.T).flat

    m3 = np.matrix([[n,                  np.sum(np.exp(  c*x))],
                    [np.sum(np.exp(c*x)),np.sum(np.exp(2*c*x))]])

    m4 = np.matrix([np.sum(y), np.sum(y*np.exp(c*x).T)])

    [a, b] = (m3.I * m4.T).flat

    return [a,b,c]