Программу, описанную в этом тексте можно скачать в разделе программы.
Вспомним, что такое аппроксимация. Допустим, у нас есть набор точек такой
как на первом рисунке:
Допустим, мы хотим провести через эти точки некоторую кривую, как на втором рисунке. Это может быть кривая вида y = A1 * x2, или y = A1 + A2*x + A3*x2 + ... Или даже y = A1 * sin(A2*x). Какой именно тип кривой надо выбрать, это дело личное. Его можно выбирать из соображений эстетики, личных предпочтений, теоретических соображений. Так вот, когда тип кривой выбран, нужно найти ряд параметров этой кривой, те самые коэффициенты A1, A2, A3...
Естественно, эти параметры должны быть такими, чтобы кривая проходила как можно ближе ко всем исходным точкам, не обязательно проходя через все точки.
Последнее требование является существенным, если его усилить и потребовать чтобы набор параметров A1, A2, A3... был таким, что для него кривая ближе всего проходит к точкам, чем для любого другого набора параметров.
Описанное выше можно формализовать так:
Пусть x1, x2, x3, ... xN - набор х-координат точек, пусть y1, y2, y3, ... yN набор у-координат тех же точек. Пусть y = f(x, A1, A2, A3, ... AK) - функция, от x, в которой есть K параметров. Рассмотрим значения функции в точках xi. Это будут числа типа f(xi, A1, A2, A3, ... AK). Понятно, что это число должно быть близко к числу yi для любого i. Значит разность
yi - f(xi, A1, A2, A3, ... AK)
должна быть малым числом, близким к нулю для любого i. И квадрат этой разности малое число.
(yi - f(xi, A1, A2, A3, ... AK))2
И сумма этих квадратов по всем i, т.е. сумма по всем точкам xi, тоже малое число
Ф = Σ(yi - f(xi, A1, A2, A3, ... AK))2
Набор A1, A2, A3, ... AK является наилучшим, если Ф - минимальна.
Формально можно составить систему уравнений, которая следует из факта минимальности Ф. Раз Ф минимальна, то частная производная по любому из параметров должна равняться нулю:
dФ/d(A1) = 0
dФ/d(A2) = 0
...
dФ/d(AK) = 0
Это ровно К уравнений для К неизвестных.
В уже описанных мною методах аппроксимации использовались очень простые функции. Это были полиномы. А при использовании полиномов в аппроксимации, для нахождения коэффициентов получается система линейных уравнений первой степени. Эта задача, хоть и не тривиальна, но довольно проста. Однако, не всегда полиномы подходят под требуемую задачу.
И тогда процесс взятия частных производных становится слишком трудоемким, а процесс решения системы уравнений чуть ли не невозможным.
В таком случае приходится использовать другие методы. Опишем один из них, который я назову "метод последовательного спуска".
Значит, мы должны найти минимум некоторой функции Ф, зависящей от К параметров. Схематически график этой функции (для Ф, зависящей от двух параметров) можно изобразить так:
Как можно найти этот минимум?
Есть метод последовательного спуска. Пусть у нас заданы некоторые начальные значения параметров А1, А2, ... АК. Их может задать человек из каких-то собственных соображений. Эти начальные параметры, к сожалению не могут быть заданы автоматически так, чтобы обеспечить нахождение максимума в любом случае. Изобразим заданные параметры на рисунке точкой:
Из этой точки нам надо как-то попасть в самый минимум функции Ф, изменяя значения параметров А1, А2,... АК. Как это сделать? Можно выбрать какое-то небольшое значение сдвига dA1, dA2,... dAK от исходной точки в направлении осей А1, А2,... АК. И посмотреть, как надо изменить параметр А1, чтобы значение функции Ф уменшилось. Т.е. надо сравнить два числа: Ф(А1+dA1, A2, A3..., AK) и Ф(А1-dA1, A2, A3..., AK). Если первое меньше второго, то надо параметр А1 увеличить до значения А1+dA1, в противном случае надо его уменьшить до значения А1-dA1. После этого надо изменять параметр А2, т.е. сравнивать числа Ф(А1, A2+dA2, A3..., AK) и Ф(А1, A2-dA2, A3..., AK). И, аналогично, изменять значения параметра А2. И так для всех параметров, вплоть до АК.
После этого, повторить процедуру.
А что делать, если например Ф(А1+dA1, A2, A3..., AK) и Ф(А1-dA1, A2, A3..., AK) больше, чем Ф(А1, A2, A3..., AK)? Ведь в этом случае после сделанного шага, значение функции Ф вырастет и мы удалимся от минимума. Это возможно, если шаг достаточно велик и мы находимся недалеко от минимума и как бы "перешагиваем" его. Решение этой проблемы существует. Для этого, если обе они больше Ф(А1, A2, A3..., AK), надо уменьшить шаг dA1, например, в два раза. И снова вычислить значения функций Ф(А1+dA1, A2, A3..., AK) и Ф(А1-dA1, A2, A3..., AK). Если они все еще больше, уменьшить шаг еще в два раза.
Когда же стоит прервать эту процедуру? Естественно это сделать тогда, когда величины шага для всех переменных станут достаточно малыми. Можно задать в программе некое очень малое число типа eps = 10-6. И когда шаг станет меньше значения параметра в eps раз, движение можно прекратить. Т.е. надо прервать движение когда dA1/A1 < eps, dA2/A2 < eps,... dAK/AK < eps.
Применим этот метод.
Допустим, в результате эксперимента у меня были получены такие значения х, у, как на следующей картинке:
Увидеть их в текстовом виде вы можете в файле spusk_data.txt. Я решил, что для аппроксимации хорошо подходит гауссообразная функция вида:
|
Пишем программу. Создадим подпрограмму, которая посчитает количество строк в файле данных. На вход этой подпрограммы надо подавать файловый указатель, который связан с этим файлом во время открытия.
int count_num_lines(FILE* InFile){
//count number of lines in input file
int nelf=0; //non empty line flag
int N = 0;
do{
nelf = 0;
while(fgetc(InFile)!='\n' && !feof(InFile)) nelf=1;
if(nelf) N++;
}while(!feof(InFile));
return(N);
}
Создадим массивы (подпрограмма allocmatrix) для хранения данных и считаем в них значения х, у (подпрограмма readmatrix). После завершения работы программы должна обязательно вызываться подпрограмма для освобождения занятой памяти (freematrix).
void allocmatrix(int N){
//allocate memory for matrixes
x = new double[N];
y = new double[N];
if(x==NULL || y==NULL){
printf("\nNot enough memory to allocate. N=%d\n", N);
exit(-1);
}
}
void freematrix(){
//allocate memory for matrixes
delete [] x;
delete [] y;
}
void readmatrix(FILE* InFile, int N){
int i=0;
//read x, y matrixes from input file
for(i=0; i<N; i++){
fscanf(InFile, "%lf", &x[i]);
fscanf(InFile, "%lf", &y[i]);
}
}
Опишем функцию Ф, которая в данной программе обозначена как f:
double f(double* param, int N){
double sum = 0;
//A - param[0]
//x0 - param[1]
//m - param[2]
for(int i=0; i<N; i++){
sum += pow(y[i] - param[0] *
exp(-(x[i]-param[1])*(x[i]-param[1])/param[2]), 2);
}
return(sum);
}
Как видно из ее описания, функция вычисляет значение функционала по значениям параметров, используя массивы значений х, у.
Теперь дело за главной функцией, где и реализован алгоритм, описанный выше:
void main(){
int Nparam = 3;
double *param, *param2;
double sum_dparam=0;
param = new double[Nparam];
param2 = new double[Nparam];
//A = 10, x0 = 5, m = 10
param[0] = 10;
param[1] = 5;
param[2] = 10;
double* dparam;
dparam = new double[Nparam];
dparam[0] = 1; //dA
dparam[1] = 1; //dx0
dparam[2] = 1; //dm
double eps = 1e-8;
double yleft, yright, y;
char InFileName[256];
FILE* InFile = NULL;
int N=0;
printf("Input file name: ");
scanf("%s", InFileName);
InFile = fopen(InFileName, "rt");
if(InFile==NULL){
printf("File not found!\n");
return;
}
N = count_num_lines(InFile);
allocmatrix(N);
rewind(InFile);
readmatrix(InFile, N);
fclose(InFile);
do{
for(int p=0; p < Nparam; p++){
y = f(param, N);
for(int pi=0; pi < Nparam; pi++){
param2[pi] = param[pi];
}
param2[p] = param[p] - dparam[p];
yleft = f(param2, N);
param2[p] = param[p] + dparam[p];
yright = f(param2, N);
if(yleft < y && yright < y){
do{
dparam[p] /= 2.0;
param2[p] = param[p] - dparam[p];
yleft = f(param2, N);
param2[p] = param[p] + dparam[p];
yright = f(param2, N);
}while(yleft > y && yright > y && dparam[p]/param[p] > eps);
}
if( yleft < y && yleft < yright){
double prev_yleft = yleft;
do{
dparam[p] *= 2;
param2[p] = param[p] - dparam[p];
yleft = f(param2, N);
}while(yleft < prev_yleft);
dparam[p] /= 2;
param[p] -= dparam[p];
}
if( yright < y && yright < yleft){
double prev_yright = yright;
do{
dparam[p] *= 2;
param2[p] = param[p] + dparam[p];
yright = f(param2, N);
}while(yright < prev_yright);
dparam[p] /= 2;
param[p] += dparam[p];
}
}
sum_dparam=0;
for(p=0; p < Nparam; p++){
sum_dparam += fabs(dparam[p]/param[p]);
}
}while(sum_dparam > eps*Nparam);
printf("\n");
for(int p=0; p < Nparam; p++){
printf("A[%d]=%lf\n", p, param[p]);
}
freematrix();
delete [] param;
delete [] param2;
}
Здесь все так, как описано в алгоритме выше, с небольшими изменениями, которые позволили ускорить работу программы. Дело в том, что часто при спуске параметры dA уменьшаются на некотором шаге, однако на следующем шаге может оказаться так что вдвое увеличенное значение dA еще не уведет нас от минимума. Поэтому, перед тем как сделать шаг в одном из выбранных направлений, программа увеличивает значение текущего параметра dA вдвое и смотрит насколько уменьшилось значение функции Ф. Если его значение увеличилось, а не уменьшилось, то это удвоение отменяется и делается шаг единичной длины.
И, наконец, надо не забыть, что для успешной компиляции программы нам требуются некоторые заголовочные файлы и определения глобальных переменных, которые мы должны написать в самом начале программы:
#include <stdio.h> #include <math.h> #include <stdlib.h> double *x=NULL, *y=NULL;
При запуске программа спрашивает, из какого файла читать данные. Надо ей написать имя файла и нажать Enter. Начальные значения параметров аппроксимации написаны в теле программы так:
param[0] = 10; //A param[1] = 5; //x0 param[2] = 10; //m dparam[0] = 1; //dA dparam[1] = 1; //dx0 dparam[2] = 1; //dm
Если вам потребуется их изменить - делайте это без страха. Я лишь рекомендую не делать стартовое значение шагов делать больше, чем 10% от значения самих параметров, иначе есть шанс очень быстро уйти далеко от зоны минимума на первых же шагах. Вы можете, конечно, изменять и количество параметров, которые задаются переменной
int Nparam = 3;и аппроксимирующую функцию f.
Покажу результаты работы программы:
A[0]=13.481637 A[1]=1.975703 A[2]=32.025858Это и есть последовательно значения параметров А, x0, m. Для доказательства того, что эти значения хорошо подходят под наш набор исходных точек, покажем график функции f(x), с этими значениями параметров и исходный набор точек:
Правда, хорошо подходят друг под друга?
