В задачах интерполяции функций по заданным значениям функции для заданного набора аргументов широко применяется формула аппроксимации функции полиномом, совпадающего в заданных точках со значениями исследуемой функции.
Обобщим эту формулу на случай функции нескольких переменных
Скалярное произведение векторов будем вычислять с учётом квадратов дисперсий значений аргументов Dx
- с целью учета равнозначности каждой координаты аргумента вычисляем Dx[k] =M2[k] − M1[k]∗M1[k]
- < A; B > = SUM A[k] * B[k] / Dx[k]
Алгоритм
1. Для заданных векторов Xi рассчитываем вектор Dx квадратов дисперсий значений аргументов
2. Для вектора X находим p ближайших точек X1,…,Xp, где расстояния между точками вычисляется с учётом квадратов дисперсий значений аргументов
3. Вычисляем значение полинома P(X) = SUM Yi*П < X-Xj, Xi-Xj > / < Xi-Xj, Xi-Xj >, где скалярное произведение вычисляется с учётом квадратов дисперсий значений аргументов
Демонстрация аппроксимации для функции двух переменных
history.txt
- -1 -1 -10
- 1 1 10
- 0 0 0
- -1 1 -8
- 1 -1 -20
Построим график значений полинома на решётке с шагом 0.1 в программе gnuplot
- predict.exe -history history.txt -input input.txt -output output.txt -p 5
- splot «output.txt»
// predict.cpp: определяет точку входа для консольного приложения. // #include "stdafx.h" struct t_previous_result { std::vector<double> x; double y; }; ///////////////////////////////////////////////////////// // Вычисление квадрата растояния между двумя векторами координат double delta(std::vector<double>& a, std::vector<double>& b, std::vector<double>& dx) { auto s = 0.0; auto i = 0; for (; i < a.size() && i < b.size() && i < dx.size(); i++) if (dx[i] > 0.0) s += (a[i] - b[i]) * (a[i] - b[i]) / dx[i]; for (; i < a.size() && i < dx.size(); i++) if (dx[i] > 0.0) s += a[i] * a[i] / dx[i]; for (; i < b.size() && i < dx.size(); i++) if (dx[i] > 0.0) s += b[i] * b[i] / dx[i]; return s; } ///////////////////////////////////////////////////////// // Вычисление растояния проекции двух векторов double scalar(std::vector<double>& a, std::vector<double>& b, std::vector<double>& dx) { auto s = 0.0; auto i = 0; for (; i < a.size() && i < b.size() && i < dx.size(); i++) if (dx[i] > 0.0) s += (a[i] * b[i]) / dx[i]; return s; } ///////////////////////////////////////////////////////// // Возвращает предсказание для указанных параметров исходя из исторических данных double predict(std::vector<double>& x, std::vector<t_previous_result>& previous_results, std::vector<double>& dx, int p) { std::vector<std::pair<t_previous_result, double>> neighbors; for (auto it = previous_results.begin(); it != previous_results.end(); ++it) { std::vector<double>& x2 = it->x; auto d = delta(x, x2, dx); std::pair<t_previous_result, double> pair(*it, d); neighbors.push_back(pair); } std::sort(neighbors.begin(), neighbors.end(), [](std::pair<t_previous_result, double> const& a, std::pair<t_previous_result, double> const& b) { return (a.second < b.second); }); neighbors.resize(std::min(p, static_cast<int>(neighbors.size()))); auto y = 0.0; for (auto iti = neighbors.begin(); iti != neighbors.end(); ++iti) { auto s = iti->first.y; std::vector<double>& xi = iti->first.x; for (auto itj = neighbors.begin(); itj != neighbors.end(); ++itj) { if (iti == itj) continue; std::vector<double>& xj = itj->first.x; std::vector<double> xxj; std::vector<double> xixj; for (auto i = 0; i < x.size() && i < xj.size(); i++) xxj.push_back(x[i] - xj[i]); for (auto i = 0; i < xi.size() && i < xj.size(); i++) xixj.push_back(xi[i] - xj[i]); s *= scalar(xxj, xixj, dx) / scalar(xixj, xixj, dx); } y += s; } return y; } ///////////////////////////////////////////////////////// // Дефолтные значения static const int _p = 3; int main(int argc, char* argv[]) { std::vector<t_previous_result> previous_results; std::vector<double> dx; char* input_file_name = nullptr; char* output_file_name = nullptr; char* previous_results_file_name = nullptr; auto p = _p; std::string line; // Поддержка кириллицы в консоли Windows // Функция setlocale() имеет два параметра, первый параметр - тип категории локали, в нашем случае LC_TYPE - набор символов, второй параметр — значение локали. // Вместо второго аргумента можно писать "Russian", или оставлять пустые двойные кавычки, тогда набор символов будет такой же как и в ОС. setlocale(LC_ALL, ""); for (auto i = 1; i < argc; i++) { if (strcmp(argv[i], "-help") == 0) { std::cout << "Usage :\t" << argv[0] << " [-input <inputfile>] [-output <outputfile>] [...]" << std::endl; } else if (strcmp(argv[i], "-input") == 0) input_file_name = argv[++i]; else if (strcmp(argv[i], "-output") == 0) output_file_name = argv[++i]; else if (strcmp(argv[i], "-history") == 0) previous_results_file_name = argv[++i]; else if (strcmp(argv[i], "-p") == 0) p = atoi(argv[++i]); } if (input_file_name != nullptr) freopen(input_file_name, "r", stdin); if (output_file_name != nullptr) freopen(output_file_name, "w", stdout); if (previous_results_file_name != nullptr) { std::vector<double> m1; std::vector<double> m2; std::ifstream history(previous_results_file_name); if (!history.is_open()) throw "Error opening file"; while (std::getline(history, line)) { std::stringstream ss(line); std::vector<double> x; std::copy(std::istream_iterator<double>(ss), std::istream_iterator<double>(), std::back_inserter(x)); auto y = x.back(); x.pop_back(); t_previous_result previous_result; previous_result.x = x; previous_result.y = y; previous_results.push_back(previous_result); for (auto i = 0; i < x.size(); i++) { if (i >= m1.size()) m1.push_back(0.0); if (i >= m2.size()) m2.push_back(0.0); m1[i] += x[i]; m2[i] += x[i] * x[i]; } } for (auto it = m1.begin(); it != m1.end(); ++it) *it /= previous_results.size(); for (auto it = m2.begin(); it != m2.end(); ++it) *it /= previous_results.size(); for (auto i = 0; i < m1.size() && i < m2.size(); i++) dx.push_back(m2[i] - m1[i] * m1[i]); } while (std::getline(std::cin, line)) { double y; std::vector<double> x; std::stringstream ss(line); std::copy(std::istream_iterator<double>(ss), std::istream_iterator<double>(), std::back_inserter(x)); y = predict(x, previous_results, dx, p); for (auto it = x.begin(); it != x.end(); ++it) std::cout << *it << ' '; std::cout << y << std::endl; } return 0; } Интернет-адрес проекта: https://github.com/dprotopopov/polylib
Используемое программное обеспечение
- Microsoft Visual Studio 2013 — среда программирования
- gnuplot — кросс-платформенный инструмент для построения графиков www.gnuplot.info
ссылка на оригинал статьи https://habrahabr.ru/post/277541/
Добавить комментарий