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

int а[4] = {2, 20, 150, 700}, Ь[4], *iEnd;

iEnd = adjacent difference(а, a+4, b, divides<int>());

МЫ получим b = {2, 15, 5, 4}, так как 30/2 = 15, 150/30 = 5 и 700/150 = 4.

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

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

118.6 118.5 118.4 118.3 118.2 118.1 118.0

-Линия регрессии у = а + Ьх

20.4 20.5 20.6 20.7 20.8 20.9 21.0 21.1 х

Рисунок 7.1. Восемь точек и вычисленная линия регрессии

Количество точек произвольно, но не может быть меньше 2. Расположение точек должно быть таким, чтобы получающаяся линия не оказалась вертикальной; в частности недопустимо, если все точки имеют одну и ту же координату х.

Координаты п точек будут находиться во входном файле, имя которого вводится пользователем. Например, восемь точек, изображенных в виде крестиков на рисунке 7.1, находятся в следующем файле data.txt:

20.5 118.1

20.6 118.25 20.65 118.2

прикладная программа; метод наименьших квадратов 207

20.75 118.4

20.8 118.4 20.85 118.5

20.9 118.45 21.0 118.5

Получающаяся линия, называемая линией регрессии, будет представлена уравнением

у = а + Ьх

где коэффициенты avib вычислены таким образом, что сумма - У(г)У (где y(J) = а + Ьх)

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

an + i X X = S у, aLx. + ЬЪх = Ё л- у

Программирование решения этой задачи с помощью STL представляется удобным, потому что:

алгоритм accumulate позволит просто вычислить значения X х и S у;,

алгоритм inner product позволит просто вычислить X и X х. у:,

последовательные контейнеры позволяют хранить переменное количество значений хиу.

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

г/i - y(i)

где, как уже упоминалось, y(i) = а + Ьх.. Если бы мы не сохраняли все значения, нам пришлось бы дважды читать входной файл. Используя тип vector<double> для хранения значений л: и г/, получаем следующую программу:

leastsq.cpp: Метод наименьших квадратов, ♦include <iostream>

ttinclude <fstream>

ttinclude <iomanip>

ttinclude <vector>

ttinclude <numeric>

ttinclude <string>

using namespace std;

typedef vector<double> array;

int ReadPoints(ifstream &file, array &x, array &y) { double xl, yl;

while (file >> xl yl)

{ x.push baclc(xl); y.push baclc(yl);

return X.size();

int CpmputeCoeff(const array &x, const array &y,

double &a, double &b) { double sx=0, sx2=0, sy=0, sxy=0;

sx = accumulate (x.begin 0 , x.endO, sx) ;

sy = accumulate (y. begin 0 , y.endO, sy) ;

sx2 = inner product (x.beginO , x.endO, X.begin0, sx2);

sxy = inner product (x.begin 0 , x.endO, y.beginO, sxy);

int n = X.size();

double D = n * sx2 - sx * sx;

if (D != 0)

{ a = (sy * sx2 - sx * sxy)/D;

b = (n * sxy - sx * sy)/D;

return 1; } else return 0;

void ShowComputedPoints(const array &x, const array &y,

double a, double b) { cout \n X yGiven yComputed\n\n

setiosflags(ios::fixed); int n = X.size(); for (int i=0; i<n; i++)

cout setw(lO) setprecision(2) x[i] setw(lO) y[i]

setw(lO) a + b * x[i] endl;

int mainO

{ array X, y;