Блог им. kurd
Если кого вдохновило сообщение smart-lab.ru/blog/680086.php, тому не обойтись без книги «NUMERICAL RECIPES. The Art of Scientific Computing. Third Edition». Качайте, пока дают
www.e-maxx-ru.1gb.ru/bookz/files/numerical_recipes.pdf
Бесплатные исходники к ней github.com/blackstonep/Numerical-Recipes
Программа svd.h из этого набора решает задачу наименьших квадратов для построения индикатора полиномиальной регрессии вместо примитивных скользящих средних.
Хорошее объяснение математической подоплёки в книге «Машинные методы математических вычислений. Форсайт, Малькольм, Моулер» en.booksee.org/book/445129
Ещё лучше — «Линейная алгебра и её применения» Гилберт Стренг
fileskachat.com/download/20151_887581203f10b39b3d7f6b84caf48a63.html
«Linear Algebra and Its Applications 4ed»
www.astronomia.edu.uy/progs/algebra/Strang- Linear_algebra_and_its_applications.pdf
Для использования программы svd.h из «NUMERICAL RECIPES» нужны тривиальные дополнения — транспонирование и перемножение матриц. Набор программ можно дополнить самодельным файлом utils.h и разместить в нём такой код:
#include <assert.h>
template <class T>
class NRdiagonal: public NRvector<T> { using NRvector<T>::NRvector; };
template <typename T>
void Multiply (const NRdiagonal<T>& a, const NRvector<T>& b
,NRvector<T>& c) {
int m = a.size();
assert (m == b.size());
c.resize (m);
for (int i = 0; i < m; ++i)
c[i] = a[i] * b[i];
}
template <typename T>
void Multiply (const NRmatrix<T>& a, const NRvector<T>& b
,NRvector<T>& c) {
int m = a.nrows(); int n = a.ncols();
assert (n == b.size());
c.resize (m);
for (int i = 0; i < m; ++i) {
c[i] = 0;
for (int j = 0; j < n; ++j)
c[i] += a[i][j] * b[j];
}
}
template <typename T>
void Transpose (const NRmatrix<T>& a, NRmatrix<T>& b) {
int m = a.nrows(); int n = a.ncols();
b.resize (n, m);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
b[i][j] = a[j][i];
}
template <typename T>
void PrintVector (char* hdr, const NRvector<T>& vec) {
cout << hdr << '\n';
for (int i = 0; i < vec.size(); ++i)
cout << " " << vec[i];
cout << '\n';
}
Пример по данным из книги Форсайта выглядит так:
#include «nr3.h»
#include «svd.h»
#include «utils.h»int main() { // Население США в млн по годам
double arg[] = { 1900, 1910, 1920, 1930, 1940, 1950, 1960, 1970 };
double val[] = { 75.99, 91.97, 105.71, 123.20, 131.67, 150.70, 179.32, 203.21 };
const int m = _countof(arg), n = 3; // макро Visual C++
double a[m][n];
for (int i = 0; i < m; ++i) {
double t = (arg[i] — arg[0]) / (arg[m-1] — arg[0]);
a[i][0] = 1;
a[i][1] = t;
a[i][2] = t * t;
}
MatDoub_IO A (m, n, &a[0][0]);
SVD svd (A);
NRdiagonal<double> d (svd.w.size());
for (int i = 0; i < svd.w.size(); ++i)
d[i] = 1 / svd.w[i];
MatDoub_O Ut;
Transpose (svd.u, Ut);
VecDoub_IO b (m, &val[0]);
VecDoub_O c (m);
Multiply (Ut, b, c);
VecDoub_O e (m);
Multiply (d, c, e);
VecDoub_O x (m);
Multiply (svd.v, e, x);
PrintVector («x», x);
Multiply (A, x, b);
puts («b»);
for (int i = 0; i < m; ++i)
printf ("%9.6f;%9.4f\n", (arg[i] — arg[0]) / (arg[m-1] — arg[0]), b[i]);
// 79.0287 90.0335 103.206 118.547 136.056 155.732 177.577 201.59
puts («b50»);
for (int i = 0; i < 50; ++i) {
double t = i / 49.0;
printf ("%9.6f;%9.4f\n", t, x[0] + t*x[1] + t*t*x[2]);
}
printf («OK»);
}
Ещё будет уместно переделать в nr3.h конструкторы, operator=, методы resize и assign для NRvector и NRmatrix так, чтобы не менять выделенную память v, если новый её размер не больше старого.
Если кто захочет содрать код, замените угловые кавычки обычными ". Углы без спроса загибает смарт-лаб.
В крайней ситуации, когда базисные функции линейно зависимы, можно показать, что матрица Р вырождена, и число обусловленности может рассматриваться как бесконечное. Следовательно, методы, избегающие высокого числа обусловленности, связанного с нормальными уравнениями, есть в то же время методы, способные обнаруживать линейную зависимость среди базисных функций. Гауссово исключение и его варианты не слишком хорошо приспособлены к выявлению такой зависимости. Наша оценка числа обусловленности в какой-то мере помогает, однако она может лишь сигнализировать об опасности, но не предложить лекарство.
Наиболее надежный метод для вычисления коэффициентов в общей задаче наименьших квадратов основан на матричной факторизации, называемой сингулярным разложением. Есть другие методы, требующие меньше машинного времени и памяти, однако они менее эффективны в том, что касается учета ошибок заданной информации, ошибок округления и линейной зависимости.»
Если Вам нужно оценить десяток к-тов, а ценовой ряд имеет длину в сотни раз больше, задача хорошо обусловлена и нет вопросов. Если Вы умудритесь, на малом числе данных вычислить сравнимое или меньшее число к-тов, то без регуляризации Вам не обойтись. Потому что SVD подход не решит проблему регуляризации сам по себе. И тем боле так, как нужно для данной задачи.
Не стал врубать.)
import matplotlib.pyplot as plt
import numpy as np
x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
z = np.polyfit(x, y, 3)
plt.plot(x,y)
plt.plot(x,z)
plt.grid()
plt.show()