#алгоритм #cpp #численные_методы
Метод наименьших квадратов
Задача написать алгоритм для нахождения коэффициентов полинома N степени с помощью
МНК по M точек. То есть даны 10 точек, и сказано аппроксимировать для полинома 3-й
степени:
y=ax^3 + a1x^2 + a2x+a3
Для этого полинома надо найти a, a1, a2, a3; если для 1 и 2 степени можно написать
вручную, то даже для 3-й степени придётся решать систему 4-х уравнений с 4 неизвестными.
Вот тут можно посмотреть, как выглядит аппроксимация для разных функций с помощью мнк.
Как запрограммировать алгоритм, который будет составлять СЛОУ и решать её?
Язык c++.
Обновление
Как в составлении так и в решении. Начал решать пока что такая система получается
для полинома N степени. Ex сумма от 1 до M (кол-во точек).
a1*Ex^n + a2*Ex^n - 1 + ... + an*Ex^0=Eyx^0
a1*Ex^n + 1 + a2*Ex^n + ... + an*Ex^1=Eyx^1
a1*Ex^2n + a2*Ex^2n - 1 + ... + an*Ex^n=Eyx^n
Даже после того как перепишу все известные Ex,Ey,Exy... в массив, то я незнаю как
решать такую СЛАУ
Решение
Прошу меня простить за кривые отступы.
Функция которая находит коэф. полинома N степени
vector* mnk(vector> pts,int n,int n2)
{
//матрица хранящая неизвестные при нужных нам коэффициентах
vector> linSys(n + 1);
double tmp=0;
vector subsid(3*n + 2);//почему 3n + 2? да потому что http://mathhelpplanet.com/static.php?p=onlayn
- mnk - i - regressionniy - analiz
for (int i=0;i<2*n + 1;i ++ )//этот цикл формирует первую половину (точнее 2/3) списка
величин x^0 x^1...x^2n + 1
{
for (int j=0;j*b=new vector;
b - >resize(n + 1);
for(int i = 0; i < n + 1; i ++ )
{
linSys[i].resize(n + 1);
(*b)[i]=subsid[2*n + 1 + i];
for(int j=0;j*x=new vector;
x - >resize(n + 1);
//это решает СЛАУ и пишет результ в x
Gauss(linSys,*b,x,n + 1);
return x;
}
Функция которая решает СЛАУ
#include
#include
#include
#include
//результат записывает в x, СЛАУ представима как A*x=B
void Gauss(vector> A, vector b, vector *x, int n)
{
std::stack > swaps;
int i, j, k, t;
double kof, s;
for (i = n - 1; i > 0; -- i) {
for (t=i, j=i - 1; j >= 0; -- j) {//Поиск максимума в столбце
if (fabs(A[i][t]) < fabs(A[i][j])) {
t = j;
}
}
if (A[i][t] == 0.0) {
return;
}
if (t != i) {
for (k = n - 1; k >= 0; -- k) {
std::swap(A[k][t],A[k][i]);
}
swaps.push(std::pair(i,t));
}
for (j=i - 1; j>=0; -- j) {//
kof = A[j][i] / A[i][i];//
A[j][i] = 0.0;//
b[j] - = b[i] * kof;
for (k = i - 1; k >= 0; -- k) {
A[j][k] - = A[i][k] * kof;
}
}
}
for (i = 0; i < n; ++ i)
{
s = 0.0;
for (j = i - 1; j >= 0; -- j) {
s + = A[i][j] * (*x)[j];
}
(*x)[i] = (b[i] - s) / A[i][i];
}
while (!swaps.empty()) {
//int j=swaps.top().first;
// int y=swaps.top().second;
std::swap((*x)[swaps.top().first], (*x)[swaps.top().second]);
swaps.pop();
}
}
Ответы
Ответ 1
Задача о полиномиальной регрессии (аппроксимации табличной функции) требует большой аккуратности, поскольку при обычных методах решения расчёт главного определителя СЛАУ (определителя Вандермонда) приводит к вычитанию близких чисел и в итоге снижает относительную погрешность вычислений, что проявляется уже для полиномов 8-10 порядка. От этих недостатков свободны методы разложения в ряды по ортогональным функциям (например, разложение в ряды Фурье по косинусам), и хотелось бы так же просто обращаться с полиномами. Однако классические ортогональные полиномы (Лежандра, Чебышёва и пр.) в дискретном варианте ортогональны только на специально подобранных неравномерных сетках. Удачный выход из такой ситуации даёт малоизвестный метод ортогональной полиномиальной регрессии (Orthogonal Polynomial Curve Fitting, Jeff Reid), в котором сначала по абсциссам заданных точек конструируются ортогональные полиномы, а затем с помощью МНК рассчитываются коэффициенты при этих полиномах. 1. ТЕОРИЯ 1.1. Постановка задачи Даны n точек { (xi,yi), i = 0, 1 ... n-1 }. Найти полином g(x) = a0 + a1x + ... + amxm с минимальной суммой ∑i = 0, 1 ... n-1 (yi - g(xi))2 среди всеx возможных полиномов n-го порядка. Решение будем искать в виде g(x) = b0 p0(x) + b1 p1(x) + ... + bm pm(x), где pj(x), j = 0, 1, ... m - семейство ортогональных полиномов, для которых должны выполняться рекуррентные соотношения pj+1(x) = (x-Aj+1) pj(x) - Bj pj-1(x), p0(x) = 1, p-1(x) = 0, (1) и условия ортогональности ∑i = 0, 1 ... n-1 pj(xi) pk(xi) = 0 при j ≠ k. (2) 1.2. Построение ортогональных полиномов. Для j = 0 формулы (1) и (2) дают: p1(x) = x - A1, ∑ (xi - A1) = 0, откуда A1 = ∑xi / n. Если для j = 0 ... k ортогональные полиномы pj уже известны, то следующий полином pk+1 должен быть ортогонален каждому из них: ∑ pj(xi) pk+1(xi) = 0. Используя рекурррентные соотношения (1), получаем: ∑ xi pj (xi) pk (xi) - Ak+1∑ pj (xi) pk (xi) - Bk∑ pj(xi) pk-1 (xi) = 0, (3) xpj(x) = pj+1(x) + Aj+1pj(x) + Bj pj-1(x), (1') Из условий (2) и (1') следует, что для j = 0 ... k-2 все три слагаемых в соотношениях (3) нулевые. При j = k - 1 условия ортогональности уничтожают второе слагаемое в формуле (3) и два младших слагаемых в (1'), поэтому Bk = ∑ p2k(xi) / ∑ p2k-1(xi). (4) В случае j = k условия ортогональности обнуляют третье слагаемое в формуле (3), поэтому Ak+1 = ∑ xip2k(xi) / ∑ p2k(xi). (5) 1.3. Построение полинома регрессии Метод наименьших квадратов для полинома g(x) = b0p0(x) + b1p1(x) + ... + bmpm(x) предполагает вычисление коэффициентов bj, j = 0 ... m, по методу наименьших квадратов, т.е. минимизацию функции невязки F(b0, b1 ... bm) = ∑ (yi - b0p0(xi) - b1p1(xi) - ... - bmpm(xi))2 по этим коэффициентам, рассматриваемым в качестве переменных. При этом в точке минимума должны выполняться необходимые условия экстремума, т.е. равенство нулю частных производных F'bk(b0, b1 ... bm): -2∑ (yi - b0p0(xi) - b1p1(xi) - ... - bmpm(xi)) pk = 0, k = 0 ... m. Разбивая каждое из этих уравнений на слагаемые и используя условия ортогональности (2), получим: ∑ yipk(xi) - bk ∑ pk2(xi) = 0, k = 0 ... m, bk = ∑ yipk(xi) / ∑ pk2(xi), k = 0 ... m. (6) Приводя полином g(x) к виду g(x) = a0 + a1x + ... + amxm, получим: aj = ∑k=j...m bk pkj, j = 0 ... m. (7) Формулы (6-7) позволяют построить алгоритм расчёта коэффициентов aj полинома g(x), если ортогональные полиномы известны. 2. АЛГОРИТМ 2.1. Значения ортогональных полиномов в узлах сетки Значения полиномов Pji = pi(xi) в узлах сетки можно вычислить по формулам (1,4-5): P0i = 1, P1i = xi - Q0 / S0, Pj+1,i = (xi - Qj / Sj) Pj,i - (Sj / Sj-1) Pj-1,i для j > 0, где Sj = ∑i Pj,i2, Qj = ∑i xi Pj,i2. Для полинома третьего порядка (m=3) алгоритм имеет вид: P0i = 1; S0 = m+1; Q0 = ∑ xi, P1i = xi - Q0 / S0; S1 = ∑ (P1i)2; Q1 = ∑ xi (P1i)2, P2i = (xi - Q1 / S1) P1i - S1 / S0; S2 = ∑ (P2i)2; Q2 = ∑ xi (P2i)2, P3i = (xi - Q2 / S2) P2i - S2i / S1i; S3 = ∑ (P3i)2; Q3 = ∑ xi (P3i)2. 2.2. Коэффициенты ортогональных полиномов Рекуррентные формулы (1,4,5) также позволяют вычислить все коэффициенты сkj ортогональных полиномов вида pk(x) = сk0 + сk1 x + ... + сkj xj +... + сkm xm, которые в общем случае представляют собой алгебраическую сумму трёх величин: сk+1,j = сk,j-1 - (Qk / Sk) сk,j - (Sk / Sk-1) сk-1,j для k = 0, 1 ... m-1, j = 0, 1, ... , k, сk,k = 1 для k = 0, 1 ... m. Первое слагаемое не следует учитывать при j = 0, третье - при j = k. Для полинома третьего порядка (m=3): с00 = 1, с10 = - Q0 / S0; с11 = 1, с20 = - (Q2 / S2) с10 - (S1 / S0) с00; с21 = с10 - (Q2 / S2) с11 - (S1 / S0) с01; с22 = 1, с30 = - (Q3 / S3) с20 - (S2 / S1) с10; с31 = с20 - (Q3 / S3) с21 - (S2 / S1) с11; с32 = с21 - (Q3 / S3) с22 - (S2 / S1) с12; с33 = 1. 2.3. Коэффициенты полинома регрессии. Коэффициенты bk в разложении g(x) = b0 p0(x) + b1 p1(x) + ... + bm pm(x) вычисляются в соответствии с (6) по формулам bk = ∑i yi Pki / Sk, k = 0 ... m. Для полинома третьего порядка: b0 = ∑i yi P0i / S0, b1 = ∑i yi P1i / S1, b2 = ∑i yi P2i / S2, b3 = ∑i yi P3i / S3. Коэффициенты aj в разложении g(x) = a0 + a1 x + ... + am xm вычисляются по формуле aj = ∑ k=j...m bk сkj, j = 0 ... m. (7) Для полинома третьего порядка (m=3): a0 = b0 с00 + b1 с10 + b2 с20 + b3 с30, a1 = b1 с11 + b2 с21 + b3 с31, a2 = b2 с22 + b3 с32, a3 = b3. 2.4. Об алгоритме Разработанный алгоритм в общем виде применим для аппроксимации заданного массива точек полиномом любого порядка при условии достаточности исходных данных и вычислительных средств. На примере полинома третьего порядка показано, что никакие СЛАУ при этом решать не требуется. 3. ТЕСТИРОВАНИЕ АЛГОРИТМА Тестирование алгоритма проведено с использованием пакета MathCAD. Для тестирования были выбраны 10-точечные массивы x и y и полиномиальная модель третьего порядка (см. скриншот). Тестирование подтвердило правильность алгоритма. Сравнение результатов расчёта с функцией regress() пакета MathCAD показало, что при одинаковой разрядности расчётов (определяемой исходными данными) алгоритм ортогональной полиномиальной регрессии точнее. 4. ПРОГРАММА НА PHP Программа сoдержит класс Ortho, конструктор которого принимает исходные массивы $x, $y и порядок модели $m и рассчитывает регрессионный полином. Для выдачи рассчитанных величин предусмотрены следующие методы: getValues() - для массива значений ортогональных полиномов в точках x; getSums() - для сумм, рассчитанных по этому массиву; getAllCoeffs() - для массивов коэффициентов (с - массивы коэффициентов ортогональных регрессионных полиномов; b - массив для расчёта регрессии по значениям ортогональных полиномов; a - коэффициенты при степенях регрессионного полинома); getRegress() - для значений регрессионного полинома и их невязок по отношению к массиву $y. /** * Polinomial regression via orthogonal regression polynomials. * Author: Yury Negometyanov * 21.03.2017 */ class Ortho { private $x; private $y; private $m; private $values; private $sums; private $coeffs; public function __construct($xx, $yy, $mm) { $this->x = $xx; $this->y = $yy; $this->m = $mm; $this->calcValuesSums(); $this->calcAllCoeffs(); } public function subtract($ar, $subtr) { if(gettype($subtr) == "array") { return array_map(function($val1, $val2){ return $val1 - $val2;}, $ar, $subtr); } else { return array_map(function($val)use($subtr){ return $val - $subtr;}, $ar); } } public function multi($ar, $mult) { if($mult == "square") { return array_map(function($a){return $a * $a;}, $ar); } elseif(gettype($mult) == "array") { return array_map(function($a, $b){return $a * $b;}, $ar, $mult); } else { return array_map(function($val)use($mult){ return $val * $mult;}, $ar); } } private function calcValuesSums() { for ($j = -1; $j < $this->m ; $j ++) { if(!isset($v)) { $v = [array_fill(0, count($this->x), 1.0)]; $s = [(float)count($this->x)]; $q = [array_sum($this->x)]; } else { $v[] = $this->subtract($this->x, end($q)/end($s)); if($j > 0) { $v[$j+1] = $this->subtract( $this->multi(end($v), prev($v)), $this->multi(prev($v), end($s)/prev($s))); } $sq = $this->multi(end($v), "square"); $s[] = array_sum($sq); $q[] = array_sum($this->multi($this->x, $sq)); } } $this->values = $v; $this->sums = ['∑P²'=>$s, '∑xP²'=>$q]; } private function calcOrthoCoeffs() { $prev = null; $ab = array_map(function($ss, $qq)use(&$prev) { if(is_null($prev)) $bb = 0; else $bb = $ss / $prev; $prev = $ss; return [$qq / $ss, $bb]; }, reset($this->sums), next($this->sums)); foreach ($ab as $k => list($aa, $bb)) { if(!isset($с)) $с = [[1.0]]; $с[] = array_fill(0, $k+2, 1.0); for ($j=0; $j <= $k; $j++) { $с[$k+1][$j] = ($j==0 ? 0.0 : $с[$k][$j-1]) - $с[$k][$j] * $aa - (($j==$k) ? 0.0 : $с[$k-1][$j] * $bb); } } return $с; } private function calcAllCoeffs() { $s = reset($this->sums); $b = []; foreach ($this->values as $key => $v) { $b[] = array_sum($this->multi($this->y, $v))/$s[$key]; } $c = $this->calcOrthoCoeffs(); $m = count($b) - 1; $a = array_fill(0, $m+1, 0.0); foreach ($a as $j => &$aj) { for ($k = $j; $k <= $m; $k++) { $aj += $b[$k] * $c[$k][$j]; } } $this->coeffs = ['c' => $c, 'b'=>$b, 'a'=>$a]; } private function calcRegress($opt ='a') { switch ($opt) { case 'a': $rev = array_reverse($this->coeffs['a']); $p = []; foreach($this->x as $xx) { $sum = 0.0; foreach ($rev as $value) { $sum *= $xx; $sum += $value; } $p[] = $sum; } break; case 'b': foreach ($this->values as $key => $arr) { if(!isset($p)) $p = $this->multi($arr, $this->coeffs['b'][$key]); else $p = $this->subtract($p, $this->multi($arr, - $this->coeffs['b'][$key])); } break; default: $p = []; break; } return $p; } public function getValues() { return $this->values; } public function getSums() { return $this->sums; } public function getAllCoeffs() { return $this->coeffs; } public function getRegress() { $values_a = $this->calcRegress('a'); $values_b = $this->calcRegress('b'); return [ 'regress_a'=>$values_a, 'regress_b'=>$values_b, 'discrep_a'=>array_sum($this->multi($this->subtract($this->y, $values_a), "square")), 'discrep_b'=>array_sum($this->multi($this->subtract($this->y, $values_b), "square")) ]; } } function print_m($text, $arr, $level=0){ $space = str_repeat(" ", $level++); echo "$space$text"; $flag = false; foreach($arr as $value) $flag = $flag || (gettype($value)=="array"); foreach($arr as $key => $value) { if(gettype($value) != "array") { echo $flag ? "
$key => $value" : " $value"; } else { print_m("
$space$key => [", $value, $level); echo " ]"; } } $level--; } $x = [0.00000, 3.36588, 3.63719, 0.56448, -3.02721, -3.83570, -1.11766, 2.62795, 3.95743, 1.64847]; $y = [3.95610, 74.84479, 89.44289, 6.46668, -14.53888, -34.55881, 1.70531, 43.80101, 109.12940, 18.81613]; /* $x = []; for ($i=0; $i < 1000; $i++) { $x[] = (float)5*sin($i); } $y = []; foreach ($x as $xx) { $y[] = (float)sin($xx/2+1); } $n = count($x); $ort = new Ortho($x, $y, $m = 80); $discrep = $ort->getRegress(); echo "Sample size = $n Model Order = $m"; echo "
Discrepancy via Powers is: {$discrep['discrep_a']}"; echo "
Discrepancy via Ortogonal Regression Polynomials is: {$discrep['discrep_b']}"; */ $ort = new Ortho($x, $y, $m=3); print_m("Issue Data:", ['Array_x'=>$x, 'Array_y'=>$y, 'Model Order'=>$m]); $values = $ort->getValues(); print_m("
The Orthogonal Regression Polynomials' Values:", $values); $sums = $ort->getSums(); print_m("
Calculated Sums:", $sums); $coeffs = $ort->getAllCoeffs(); print_m("
Calculation of Polynomials' Coefficients:", $coeffs); $discrep = $ort->getRegress(); print_m("
Regression via Powers (a) & via Orthogonal Regression Polynomials (b) :", $discrep); 5. РЕЗУЛЬТАТЫ. Результаты на тестовой выборке: Issue Data: Array_x => [ 0 3.36588 3.63719 0.56448 -3.02721 -3.8357 -1.11766 2.62795 3.95743 1.64847 ] Array_y => [ 3.9561 74.84479 89.44289 6.46668 -14.53888 -34.55881 1.70531 43.80101 109.1294 18.81613 ] Model Order => 3 The Orthogonal Regression Polynomials' Values: 0 => [ 1 1 1 1 1 1 1 1 1 1 ] 1 => [ -0.782083 2.583797 2.855107 -0.217603 -3.809293 -4.617783 -1.899743 1.845867 3.175347 0.866387 ] 2 => [ -7.2281734230875 2.8073623836192 4.6030924341892 -7.1264829728247 3.0992779145833 8.9585998726813 -5.5494580486592 -1.3320551385999 6.9121153510662 -5.1442783729677 ] 3 => [ 3.6796621840027 -2.8171823354639 3.1956859234941 -3.0255967835539 8.7399448097287 -12.367028020194 15.20316889478 -12.281110609376 12.297473127759 -12.625017191178 ] Calculated Sums: ∑P² => [ 10 69.17098425001 328.7768235086 962.75401849569 ] ∑xP² => [ 7.82083 -27.512890311149 -1.7129366948718 250.39738237352 ] Calculation of Polynomials' Coefficients: c => [ 0 => [ 1 ] 1 => [ -0.782083 1 ] 2 => [ -7.2281734230875 -0.38433110143359 1 ] 3 => [ 3.6796621840027 -11.983278954678 -0.37912107270775 1 ] 4 => [ 20.209167857339 7.9217601883617 -14.812965853531 -0.63920555697172 1 ] ] b => [ 29.906462 15.897655646369 2.3791287087584 1.000001267701 ] a => [ 3.9560877250835 2.9999883433859 2.0000071554385 1.000001267701 ] Regression via Powers (a) & via Orthogonal Regression Polynomials (b) : regress_a => [ 3.9560877250835 74.844667502068 89.443009253453 6.4666635861522 -14.538829417751 -34.558843534167 1.7053151756382 43.801163134951 109.12933595098 18.816050623593 ] regress_b => [ 3.9560877250835 74.844667502068 89.443009253453 6.4666635861522 -14.538829417751 -34.558843534167 1.7053151756382 43.801163134951 109.12933595098 18.816050623593 ] discrep_a => 6.7210313148693E-8 discrep_b => 6.7210313149619E-8 Результаты работы программы на тестовой выборке соответствуют ожидаемым. На выборке в n=1000 точек при порядке модели 80 получено: Sample size = 1000 Model Order = 80 Discrepancy via Powers is: 0.3261862533562 Discrepancy via Ortogonal Regression Polynomials is: 4.2188699484271E-26 Таким образом, расчёт регрессии через ортогональные полиномы более устойчив по отношению к ошибкам округления. 6. ВЫВОДЫ. Алгоритм можно рекомендовать как удачную альтернативу традиционным алгоритмам для полиномиальной регрессии произвольного порядка. Исходные данные и полученные результаты могут быть использованы для отладки на других языках программирования.
Комментариев нет:
Отправить комментарий