Аппроксимация экспонентой , нахождение коэффициентов

[Super-Jeronimo]

Всем привет. Нужна помощь.
Необходимо аппроксимировать полученные табличные данные экспонентой.
То есть функцией вида: f(x)=a*exp(b*x)+c .
Допустим, в Маткаде это делается с помощью функции expfit - находятся коэффициенты a, b и c.
Проблема в том, что мне нужен сам алгоритм нахождения этих трех коэффициентов без участия Маткада и любой другой подобной программы.
Заранее спасибо за помощь.

Добавлено 27.08.05, 13:46
Кстати, вот еще что:
где-то прочитал, что эти же коэффициенты используются в методе Левенберга-Марквардта.
Может поможет?

[esperanto]

Существует бесконечное множество способов, сделать нужную вам апроксимацию.


Например, Вы , можете воспользоваться методом наименьших квадратов.

[shadeofgray]

Цитата esperanto @ 27.08.05, 18:55

Например, Вы , можете воспользоваться методом наименьших квадратов

Только нелинейным. Линейный здесь не работает.

[Super-Jeronimo]

Хм...
А нельзя ли по-подробнее рассказать.

[shadeofgray]

Ну, есть линейный МНК, описанный например здесь - http://alglib.sources.ru/interpolation/genleastsq.php

Если бы было не exp(b*x), а exp(x), то все параметры (a, c) входили бы в выражение линейно, и можно было бы применить линейный МНК, сводящий задачу к линейным уравнениям. Здесь же параметр b входит нелинейно, так что для аппроксимации по МНК придется использовать более сложный алгоритм (Левенберга-Марквардта), решающий нелинейные уравнения.

Линейный алгоритм можно скачать по ссылке выше, но нелинейного алгоритма там нет. Вообще, лучше бы использовать специализированный под эту задачу алгоритм, но если надо по-быстрому получить решение, то можно сделать так - скачать алгоритм минимизации функции методом главных осей http://alglib.sources.ru/extremums/principalaxis.php и минимизировать функцию
G(a,b)=(f(x1)-a*exp(b*x1)-c)^2 + (f(x2)-a*exp(b*x2)-c)^2 + ...

Некрасиво, но работает.

P.S. Мне почему-то кажется, что где-то существует специализированный алгоритм для решения такой задачи, но никак не могу вспомнить детали.

[Super-Jeronimo]

Спасибо за информацию.
Есть еще вопрос: где можно взять описание алгоритма Левенберга-Марквардта на русском языке. А то нашел на английском, да там мало чего понятно.

Если уж не получится с Левенбергом-..., то придется реализовывать так, как предложили вы.

[shadeofgray]

Честно говоря, не представляю. В рунете вообще относитель мало научной информации, что неудивительно, если сравнить число ученых, говорящих по-английски, и число ученых, говорящих по-русски.

[Super-Jeronimo]

Ок, ясно. Будем искать.
Но все равно спасибо за предоставленную информацию.
Как только все получится - сообщу

[Древень]

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

[Super-Jeronimo]

Экспонента получается (по-крайней мере должна получится) убывающей, либо вогнутой, либо выпуклой, неважно. Тем более влияние коэффициента 'C' можно посмотреть потом.

To Древень:
не подскажешь, как можно реализовать то, что ты предложил.

[Super-Jeronimo]

To shadeofgray

Пытаюсь реализовать задуманное с помощью вот этого. (http://alglib.sources.ru/extremums/principalaxis.php)
Но что-то не получается.
Мои действия:
x - вектор из 3-х параметров: a, b, c.
Создаю 2 массива с исходными точками: по Х и по У.
Функцию f(x) определил след. образом:

double f(ap::real_1d_array& x)
{
    int i = 1;
    int point;
    double res = 0;
    for(i = 1; i <= M; i++)
    {
        point=*pointsY[i];
        res = res + (point-x(1)*exp(x(2)*point)-x(3));
    }
    return res;
}

В main() заполняю 2 массива с точками, а потом вызываю principalaxisminimize(....).
В итоге получаю "Segmentation fault"
Что не так?

[shadeofgray]

Во-первых, какая-то неправильная формула получается. Вы точно уверены, что в теле цикла вы вычисляете "G(a,b,c)=(f(x1)-a*exp(b*x1)-c)^2 + (f(x2)-a*exp(b*x2)-c)^2 + ..." ?

Вот это

point=*pointsY[i];
res = res + (point-x(1)*exp(x(2)*point)-x(3));

больше похоже на (f-a*exp(b*f)-c). И где возведение отклонений в квадрат? Вы даже не модули отклонений суммируете, а сами отклонения.

По поводу segmentation fault. Будет лучше, если вы приведете код, вызывающий principalaxisminimize, но даже не видя его, я подозреваю, что вы забыли вызвать метод setbounds у динамического массива X, передаваемого в функцию, и задать начальное приближение к решению. Если же нет, то будем разбираться.

Добавлено 31.08.05, 06:47
Только что посмотрел вашу задачу на компе. Сходимость к решению очень сильно зависит от выбора начального приближения, алгоритм постоянно "проваливается" в один и тот же локальный минимум. Рекомендую использовать повторные запуски алгоритма со случайных позиций с выбором наилучшего из полученных результатов.

[Super-Jeronimo]

Хм.
Согласен, протерял квадрат в формуле.
Нашел косяк в функции:
"G(a,b,c)=(f(x1)-a*exp(b*x1)-c)^2 + (f(x2)-a*exp(b*x2)-c)^2 + ..."
Теперь она считается вот так:

double f(ap::real_1d_array& x)
{
    int i = 1;
    int X_point, Y_point;
    double res = 0;
    for(i = 1; i <= M; i++)
    {
        X_point=*pointsX[i];
        Y_point=*pointsY[i];
        res = res + (Y_point-x(1)*exp(x(2)*X_point)-x(3))*(Y_point-x(1)*exp(x(2)*X_point)-x(3));
    }
    return res;
}

Добавил в main() setbounds для X.
Все равно в итоге получаю "Segmentation fault".
Пытаюсь разобраться.....

[shadeofgray]

Приведите исходный код, вызывающий principalaxisminimize. Я только что решил эту задачу у себя на компе, и всё работает нормально. И где именно возникает segmentation fault - в principalaxisminimize, до её вызова или после?

[Super-Jeronimo]

Все, с Segmentation fault разобрался.
Но значения получаются не очень удовлетворительными.
При начальных a=b=c=1 получается следующее:
A=589.950938, B=-0.075877, C=-170.350994
Minimize_Result= 2411.259707

Вот если аппроксимировать те же данные в маткаде на странице
http://twt.mpei.ac.ru/mas/worksheets/Fit_f_x_a_b_c.mcd
Х = 6 7 8 9 10 11
У = 209 154 161 161 80 0
F = a*e^(b*x)+c
a=1, b=1, c=1,
то значения отличаются значительно.

[shadeofgray]

Да, я это уже заметил. Я точно не знаю, с чем это связано. Возможно, сходимость метода Левенберга-Марквардта на тех задач, для решения которых он предназначен, значительно лучше сходимости неспециализированных методов.

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

Было бы очень интересно узнать, каким же именно алгоритмом пользуются на том сервере.

[Super-Jeronimo]

Да я вот тоже все пытался узнать. Но там у них стоит Mathcad Application Server и конкретно в этом примере (который я приводил) используются встроенные функции. Для моего случая используется как раз функция Expfit(). Единственное, что удалось "раскопать", так это то, что коэффициенты считаются по Левенбергу-.....

Но все равно спасибо. Какое-никакое решение все же есть.
Премного благодарен.

[Древень]

Цитата Super-Jeronimo @ 29.08.05, 08:15

не подскажешь, как можно реализовать то, что ты предложил.

То есть если всё-таки попробовать без минимизации:

Сначала приближённо определяем константу c:
Предполагаем, что экпонента убывает достаточно быстро. Тогда на правом краю отрезка значение функции приблизительно равно c.

log(y - c) = log(a) + b*x

То есть теперь имеем набор значений:
y' = log(y - c)
и
x

Которые нужно связать линейно, то есть найти коэффициенты

Z1 = log(a)
Z2 = b

Как это сделать? Ну для линейной аппроксимации наверняка есть точные аналитические формулы. А по простому - можно, например, вычислить приближённое значение данных y' на краях отрезка (усреднив несколько чисел вблизи краёв). А потом через эти две точки "провести" прямую.


А теперь, с учётом полученных a и b можно более точно вычислить c.

Сообщение отредактировано: Древень - 02.09.05, 05:17

[shadeofgray]

Цитата Древень @ 02.09.05, 05:15

Предполагаем, что экпонента убывает достаточно быстро

А если она возрастает? Или недостаточно быстро?

[amk]

Если возрастает, то ситуация - симметричная описанной, для определения c смотрим налево.
Если недостаточно быстро, то возникнут проблемы. Точность и так будет малой, а тут ещё прицепиться не к чему.

[Super-Jeronimo]

Всем привет.
Вот нашел еще один алогритм. Может кому пригодится:
"Решение задачи безусловной минимизации диффеpенциpуемой функции многих переменных, представимой в виде суммы квадратов, квазиньютоновским методом с использованием пpоцедуpы Левенбеpга - Маpкуаpдта."
http://www.srcc.msu.su/num_anal/lib_na/cat/mn/mnavr.htm

Только пока еще сам не разобрался как там надо задавать данные и функцию.

А вот здесь перечислены все доступные алгоритмы:
http://www.srcc.msu.su/num_anal/lib_na/int_m/int_m2.htm

[shadeofgray]

Если окажется, что у него сходимость лучше, скажете, ок?

А ещё только что подумалось, что вместо задачи трехмерной оптимизации можно решать задачу ОДНОМЕРНОЙ оптимизации. Если зафиксировать параметр b, то вместо нелинейной задачи МНК мы получаем линейную задачу размерности два, которая сводится к решению системы линейных уравнений. Ну так вот, надо проводить минимизацию по b, вычисляя параметры a и c для каждого конкретного b.

Одномерная оптимизация значительно легче многомерной, так что...

[Super-Jeronimo]

Ну вообщем так:
Более-менее разобрался с тем, как вводятся данные и функции, но видимо у меня где-то косяк, так как результаты получаются оооооочень хреновые (почти равны начальным приближениям), да и итераций происходит всего 0-2.
Так что о сходимости пока ничего не смогу сказать. Как только что-нибудь получится - сообщу.

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

[shadeofgray]

Цитата Super-Jeronimo @ 06.09.05, 02:08

но в моем случае это займет больше времени при переборе фиксированного b.

Не думаю. Одно вычисление a и c при известном b по времени близко к нескольким вычислениям значения целевой функции. Я думаю, что экономия будет.

[Super-Jeronimo]

Короче, разобрался вот с этим алгоритмом:
http://www.srcc.msu.su/num_anal/lib_na/cat/mn/mnavr.htm
Точнее не с самим алгоритмом, а со способом задания функции и данных (косяк исправил).
Главное аккуратно выставить размеры всех массивов и матриц.
Сходится просто зашибись, результат получается примерно такой же, как и в проге Маткада
( http://twt.mpei.ac.ru/mas/worksheets/Fit_f_x_a_b_c.mcd ), различия в основном начинаются с третьего порядка, но это можно исправить, поиграв с различными эпсилонами и дельтами.
Вот код, который необходим для аппроксимации данных (5 точек) функцией a*e^(b*x)+c:

#include "f2c.h"
#include "mnavr_c.c"
#include <math.h>
 
#define dimN 3
#define dimM 5
#define dim_XJTJ 6
 
int main(void)
{
    /* Initialized data */
    static int m = dimM;
    static float x[dimN] = { 1,1,1 };
    static int n = dimN;
    static int nsig = 6;
    static float eps = 1e-6f;
    static float delta = 1e-4f;
    static int maxfn = 100;
    static int iopt = 0;
    static int ixjac = dimN;
    static float parm[4] = { 0.f, 0.f, 0.f, 0.f };
 
    /* Local variables */
    static float xjac[dimM*dimN];
    extern int func_c();
    static float xjtj[dim_XJTJ], work[31], f[dimM];
    static int j, infer;
    extern int mnavr_c(U_fp, int *, int *, int *, float *, float *, int *,
                       int *, float *, float *, float *, float *, float *,
                       int *, float *, float *, int *, int *);
    static int ier;
    static float ssq;
 
    mnavr_c((U_fp)func_c, &m, &n, &nsig, &eps, &delta, &maxfn, &iopt, parm,
            x, &ssq, f, xjac, &ixjac, xjtj, work, &infer, &ier);
    printf("\nIERR=%5i", ier);
    printf("\nINFER=%5i", infer);
    printf("\nSSQ=%16.7e", ssq);
    for (j = 0; j < dimN; j++)
        printf("\n X[%d]=%f", j, x[j]);
    for (j = 0; j < dimM; j++)
        printf("\n F[%d]=%f", j, f[j]);
    for (j = 1; j <= 5; ++j) {
         printf("\nWORK[%d]=%16.7e", j, work[j-1]);
    }
    return 0;
} /* main */
 
int func_c(float *x, int *m, int *n, float *f)
{
    /* Parameter adjustments */
    --f;
//    --x;
 
    /* Function Body */
    f[1] = 209.f-x[0]*exp(6.f*x[1])-x[2];
    f[2] = 154.f-x[0]*exp(7.f*x[1])-x[2];
    f[3] = 161.f-x[0]*exp(8.f*x[1])-x[2];
    f[4] = 161.f-x[0]*exp(9.f*x[1])-x[2];
    f[5] = 80.f-x[0]*exp(10.f*x[1])-x[2];
 
    return 0;
} /* func_c */

Может быть кому-нибудь пригодится.

[shadeofgray]

Спасибо. А можно протестировать программу и сообщить результаты при параметре iopt=1? Мне интересно, как сходится именно схема Левенберга-Марквардта, а не схема Брауна?

[Super-Jeronimo]

---------------------------
При iopt=0:
a=-0.002755
b=1.051729
c=184.126968

Минимум функции "G(a,b,c)=(f(x1)-a*exp(b*x1)-c)^2 + (f(x2)-a*exp(b*x2)-c)^2 + ..." =
=1635.983765

----------------------------
При iopt=1:
a=-0.005601
b=0.981501
c=185.619095

Минимум функции "G(a,b,c)=(f(x1)-a*exp(b*x1)-c)^2 + (f(x2)-a*exp(b*x2)-c)^2 + ..." =
=1637.504761

----------------------------
Сходится лучше по Брауну при равных условиях, хотя на другом тесте может быть и наоборот. Надо смотреть