Ядерное сглаживание

Использование ядерной регрессии для прогноза спроса в сетевых магазинах

Ядерное сглаживание

Доброго времени суток, уважаемые хабровчане! В данной публикации речь пойдет о модели прогноза спроса на товары в сетевых магазинах и ее реализации на C++.

Постановка задачи

Допустим, у нас имеется сеть магазинов, в каждый из которых завозят товары. Товары (для модели прогноза) попадают в каждый магазин произвольным образом. За некий период времени мы имеем статистику — сколько в каждом магазине продано тех или иных товаров.

Требуется спрогнозировать продажи товаров за период времени, аналогичный выбранному, для всех магазинов по всем товарам, которые в них не завозились. Примечания и допущения постановки задачи

  • Товары, завезенные в магазины, не заканчивались за период сбора статистики.
  • Если в магазин завезли новые для него товары (при том, что старые товары остались), продажи не перераспределяться между старыми и новыми товарами. Статистика по старым товарам останется прежней, просто кто-то дополнительно покупает новые товары. Прогнозирование при невыполнении этого условия потребует дополнительных данных о том, как насыщается спрос при увеличении количества товаров.
  • Период, за который собирали статистику, и период, для которого нужно сделать прогноз, идентичны по спросу.

Метод решения задачи

Задача решается с помощью метода ядерной регрессии (используя формулу Надарая-Ватсона). В качестве функции ядра используется квадратическая ядерная функция с шириной окна $inline$kernel_H$inline$. «Расстояние» между товарами с разной ценой считается как:

$$display$$R_{prod}=\left|log_{10}\left(\frac{price_1}{price_2}\right)\right|$$display$$

«Расстояние» между двумя магазинами считается как произведение задаваемого вручную коэффициента $inline$K_{shops}$inline$ и отрицательного десятичного логарифма взвешенной корреляции. Взвешенная корреляция вычисляется как произведение отношения количества типов товаров, которые есть в обоих магазинах к общему количеству типов товаров (вес «правдоподобия») и корреляции количеств проданных товаров (тех, которые общие для этих магазинов):

$$display$$Weight_{shops}=\frac{N_{1,2}}{N_{total}}$$display$$

$$display$$R_{shops} = – K_{shops} \cdot log_{10} \left( Correlation_{shops} \cdot Weight_{shops} \right) $$display$$

Расстояние между двумя товарами в двух разных магазинах считается как корень из суммы квадратов расстояний между данными товарами и между данными магазинами. Задаваемый вручную коэффициент для «расстояния» между двумя магазинами совместно с задаваемой шириной окна позволяют учитывать «расстояния» между товарами и магазинами в нужной нам пропорции. Я запрограммировал два метода прогнозирования. «Общий» метод учитывает все проданные товары во всех магазинах. «Крестообразный» метод учитывает все проданные товары в текущем магазине и количества проданных прогнозируемых товаров в других магазинах.

Разница в результатах прогнозирования разными методами есть, но, на первый взгляд, небольшая (проверялась при окне $inline$kernel_H=3.0$inline$, коэффициенте расстояния магазинов $inline$K_{shops}=1.0$inline$ на матрице 20 магазинов * 20 товаров). При этом, при больших окнах функции ядра второй метод значительно быстрее первого.

Исходный код

Статистика представляет собой таблицу, в которой записаны количества проданных товаров (столбцы соответствуют товарам, строки — магазинам). Если товар не поступал в магазин, в соответствующей ячейке таблицы располагается “-1”. Для удобства в первой строке таблицы указаны цены товаров. Заголовочный файл класса хранения данных DataHolder.

h#pragma once#include class DataHolder{protected: // prices vector std::vector prices; // mask of initial data: 1 if has init and 0 if not std::vector maskData; // initial data matrix std::vector initData; // full data matrix (with prognosis) std::vector fullData; public: DataHolder() {} virtual ~DataHolder() {} // Load init data from csv file // calculate range between shops // @param filename – file path of loading data void loadData(std::string filename); // Save full data in csv file // @param filename – file path of saving data void saveData(std::string filename);};

Исходный файл класса хранения данных DataHolder.cpp#include “DataHolder.h” void DataHolder::loadData(std::string filename){ // Here is loading data from filename to initData fullData = initData;} void DataHolder::saveData(std::string filename){ // Here is saving data from fullData to filename}

Заголовочный файл класса расчета расстояний Prognoser.h#pragma once#include #include #include “DataHolder.h”class RangeCalculator : public DataHolder{ friend class Prognoser;private: // count of shops int shopsCount; // count of products int productsCount; // marix of range squares between prices std::vector priceRanges; // marix of range squares between shops std::vector shopsRanges; // sums of earn of shops std::vector shopEarnSums; public: RangeCalculator() {} virtual ~RangeCalculator() {} // calculate ranges and advanced parameters void calculateRanges();private: // calculate and fill priceRanges void fillPriceRanges(); // calculate and fill shopsRanges void fillShopsRanges(); // calculate range between shops // @param i – first shop index // @param j – second shop index double calculateShopsRange(int i, int j); // calculate earnings of shop // @param i – shop index double calculateShopEarnings(int i); // calculate shopEarnSums void calculateShopEarnSums(); // calculate range between different shops // @param i – first shop index // @param j – second shop index double calculateShopsRangeByCorrelation(int i, int j); // calculate correlation between two vectors double calculateCorrelation(const std::vector &X, const std::vector &Y);};

Исходный файл класса расчета расстояний Prognoser.cpp#include “RangeCalculator.h” void RangeCalculator::calculateRanges(){ shopsCount = initData.size(); productsCount = initData[0].size(); calculateShopEarnSums(); fillPriceRanges(); fillShopsRanges();} void RangeCalculator::fillPriceRanges(){ // initialize vector priceRanges.resize(productsCount); for (int i = 0; i < productsCount; ++i) priceRanges[i].resize(productsCount); // calculate and fill double range = 0, range2 = 0; for (int i = 0; i < productsCount; ++i) { priceRanges[i][i] = 0; for (int j = i + 1; j < productsCount; ++j) { range = log10((double)prices[i] / (double)prices[j]); range2 = range * range; priceRanges[i][j] = range2; priceRanges[j][i] = range2; } }} void RangeCalculator::fillShopsRanges(){ // initialize vector shopsRanges.resize(shopsCount); for (int i = 0; i < shopsCount; ++i) shopsRanges[i].resize(shopsCount); // calculate and fill double range = 0, range2 = 0; for (int i = 0; i < shopsCount; ++i) { shopsRanges[i][i] = 0; for (int j = i + 1; j < shopsCount; ++j) { range = calculateShopsRange(i, j); range2 = range * range; shopsRanges[i][j] = range2; shopsRanges[j][i] = range2; } }} double RangeCalculator::calculateShopsRange(int i, int j){ if (i != j) { return calculateShopsRangeByCorrelation(i, j); } else return 0;} double RangeCalculator::calculateShopsRangeByCorrelation(int i, int j){ // collects products of shops, that are in both shops std::vector maskX = maskData[i]; // mask of products in first shop std::vector maskY = maskData[j]; // mask of products in second shop std::vector maskXY(productsCount); // mask of products in both shops for (int k = 0; k < productsCount; ++k) maskXY[k] = maskX[k] * maskY[k]; // count of products in first shop int vecLen = std::accumulate(maskXY.begin(), maskXY.end(), 0); // weight of correlation calculation double weightOfCalculation = (double)vecLen / (double)productsCount; // calculating X and Y vectors std::vector X(vecLen); std::vector Y(vecLen); int p = 0; for (int k = 0; k < productsCount; ++k) { if (maskXY[k]) { X[p] = initData[i][k]; Y[p] = initData[j][k]; ++p; } } // calculating range between shops double correlation = calculateCorrelation(X, Y); // correlation double weightedCorrelation = correlation * weightOfCalculation; double range = -log10(fabs(weightedCorrelation) + 1e-10); return range;} double RangeCalculator::calculateCorrelation(const std::vector &X, const std::vector &Y){ int count = X.size(); double sumX = (double)std::accumulate(X.begin(), X.end(), 0); double sumY = (double)std::accumulate(Y.begin(), Y.end(), 0); double midX = sumX / (double)count; double midY = sumY / (double)count; double cor1 = 0, cor2 = 0, cor3 = 0; for (int i = 0; i < count; ++i) { cor1 += (X[i] - midX) * (Y[i] - midY); cor2 += (X[i] - midX) * (X[i] - midX); cor3 += (Y[i] - midY) * (Y[i] - midY); } double cor = cor1 / sqrt(cor2 * cor3); return cor;} double RangeCalculator::calculateShopEarnings(int i){ double earnings = 0; for (int j = 0; j < productsCount; ++j) { earnings += prices[j] * initData[i][j]; } return earnings;} void RangeCalculator::calculateShopEarnSums(){ shopEarnSums.resize(shopsCount); for (int i = 0; i < shopsCount; ++i) shopEarnSums[i] = calculateShopEarnings(i);} Заголовочный файл класса расчета прогнозов Prognoser.h#pragma once#include “RangeCalculator.h” class Prognoser{private: // shared_ptr by RangeCalculator object std::shared_ptr rc; // koefficient to multiply by shops ranges double kShopsR; // square of kernel window width double h2; public: // constructor // @param rc – shared_ptr by RangeCalculator object // @param kShopsR – koefficient to multiply by shops ranges // @param h – kernel window width Prognoser(std::shared_ptr rc, double kShopsR, double h); virtual ~Prognoser(); // prognose all missing data // @param func_ptr – calculating of weighted sums function pointer void prognose(void (Prognoser::*func_ptr)(int, int, double&, double&)); // calculate weighted sums with “cross” method // @param shopInd – shop index // @param prodInd – product index // @param weightsSum – sum of weights // @param contribSum – weighted sum of contributions void calculateWeightSumsCross(int shopInd, int prodInd, double& weightsSum, double &contribSum); // calculate weighted sums with “total” method // @param shopInd – shop index // @param prodInd – product index // @param weightsSum – sum of weights // @param contribSum – weighted sum of contributions void calculateWeightSumsTotal(int shopInd, int prodInd, double& weightsSum, double &contribSum);private: // calculate kernel // @param r2h2 – r*r/h/h, where r – range and h – window width double calculateKernel(double r2h2); // calculate prognosis of selected position with selected method // @param shopInd – shop index // @param prodInd – product index // @param func_ptr – calculating of weighted sums function pointer int calculatePrognosis(int shopInd, int prodInd, \ void (Prognoser::*func_ptr)(int, int, double&, double&));};

Исходный файл класса расчета расстояний Prognoser.cpp#include “Prognoser.h” Prognoser::Prognoser(std::shared_ptr rc, double kShopsR, double h){ this->rc = rc; this->kShopsR = kShopsR; this->h2 = h * h;} Prognoser::~Prognoser(){} double Prognoser::calculateKernel(double r2h2){ if (r2h2 > 1) return 0; else return (1 – r2h2) * (1 – r2h2) * 15 / 16;} void Prognoser::prognose(void (Prognoser::*func_ptr)(int, int, double&, double&)){ for (int i = 0; i < rc->shopsCount; ++i) { for (int j = 0; j < rc->productsCount; ++j) { if (!rc->maskData[i][j]) { rc->fullData[i][j] = calculatePrognosis(i, j, func_ptr); } } }} void Prognoser::calculateWeightSumsCross(int shopInd, int prodInd, double& weightsSum, double &contribSum){ double r2 = 0, r2h2 = 0, weight = 0; // calculate sums by shops for (int i = 0; i < rc->shopsCount; ++i) { if (rc->maskData[i][prodInd] && shopInd!=i) { r2 = rc->shopsRanges[shopInd][i]; r2 = r2 * kShopsR * kShopsR; r2h2 = r2 / h2; weight = 0 ? r2h2 >= 1. : calculateKernel(r2h2); weightsSum += weight; contribSum += weight * rc->initData[i][prodInd] * \ rc->shopEarnSums[i] / rc->shopEarnSums[shopInd]; } } // calculate sums by products for (int j = 0; j < rc->productsCount; ++j) { if (rc->maskData[shopInd][j] && prodInd != j) { r2 = rc->priceRanges[prodInd][j]; r2h2 = r2 / h2; weight = 0 ? r2h2 >= 1. : calculateKernel(r2h2); weightsSum += weight; contribSum += weight * rc->initData[shopInd][j]; } }} void Prognoser::calculateWeightSumsTotal(int shopInd, int prodInd, double& weightsSum, double &contribSum){ double r2 = 0, r2h2 = 0, weight = 0; for (int i = 0; i < rc->shopsCount; ++i) { for (int j = 0; j < rc->productsCount; ++j) { if (i != shopInd || j != prodInd) if(rc->maskData[i][j]) { r2 = rc->shopsRanges[shopInd][i]; r2 = r2 * kShopsR * kShopsR; r2 += rc->priceRanges[prodInd][j]; r2h2 = r2 / h2; weight = 0 ? r2h2 >= 1. : calculateKernel(r2h2); weightsSum += weight; contribSum += weight * rc->initData[i][j] * \ rc->shopEarnSums[i] / rc->shopEarnSums[shopInd];; } } }} int Prognoser::calculatePrognosis(int shopInd, int prodInd, \ void (Prognoser::*func_ptr)(int, int, double&, double&)){ double weightsSum = 0; // sum of weights double contribSum = 0; // sum of weighted contributions //calculateWeightSumsCross(shopInd, prodInd, weightsSum, contribSum); (this->*func_ptr)(shopInd, prodInd, weightsSum, contribSum); int prognosis = -1; if (weightsSum > 0) prognosis = int(contribSum / weightsSum); return prognosis;}

Ссылки

К. В. Воронцов. Лекции по алгоритмам восстановления регрессии.

  • C++
  • машинное обучение
  • big data
  • ядерная регрессия
  • алгоритмы

Хабы:

  • Занимательные задачки
  • C++
  • Big Data
  • Математика
  • Машинное обучение
  • 23 июня 2017 в 11:16
  • 20 февраля 2017 в 16:01
  • 15 февраля 2015 в 22:06

Источник: https://habr.com/ru/post/359170/

Ядерная оценка неизвестной плотности вероятности

Ядерное сглаживание

Совершенствование языка MQL5 в плане его быстродействия и постоянный рост производительности персональных компьютеров привели к тому, что пользователи платформы MetaTrader 5 все чаще для анализа рынка стали использовать достаточно сложные, развитые математические методы. Эти методы могут принадлежать различным областям экономики, эконометрики, статистики, но в любом случае при их использовании приходится сталкиваться с таким понятием как функция плотности вероятности.

Многие распространенные методы анализа разрабатывались исходя из предположения о нормальности распределения самих данных или о нормальности ошибок в используемой модели.

Кроме того, часто для оценки результатов анализа важно знать распределение различных компонент используемой модели.

И в том, и в другом случае возникает задача создания, по возможности универсального “инструмента”, позволяющего производить оценку неизвестной плотности вероятности.

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

Первоначальная идея заключалась в том, чтобы при этом не использовать никаких внешних средств, то есть предполагалось, что все должно быть реализовано только средствами MQL5. Но в результате от этой первоначальной идеи пришлось несколько отклониться.

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

Непосредственно вычисление самой оценки и ее визуализация, то есть отображение в виде графика или диаграммы. Вычисления, конечно, реализованы средствами MQL5, а визуализацию пришлось реализовать путем создания HTML-странички и отображением ее в WEB-браузере. Такое решение было выбрано в связи с желанием получить в конечном итоге графический результат в векторной форме.

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

Тем более что в недалеком будущем можно ожидать появления в составе MetaTrader 5 различного рода библиотек, наверное, в том числе и графических (насколько известно, работа в этом направлении сейчас ведется).

В случае появления в составе MetaTrader 5 развитых средств, предназначенных для построения графиков и диаграмм, переделать в предлагаемом решении отображающую часть, при желании будет несложно.

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

Хотя предложенное в статье решение и не является узкоспециализированным, но и назвать его действительно полностью универсальным нельзя.

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

Поэтому, если имеется предварительная информация о характере оцениваемого распределения, то всегда можно выбрать решение, которое будет наиболее подходящим для данного конкретного случая.

Но, тем не менее, в дальнейшем будем считать, что о характере оцениваемой плотности нам ничего не известно.

Наверняка такой подход в какой-то степени повлияет на качество но будем надеяться, что это окупится за счет возможности производить оценки совершенно различных по характеру плотностей.

В связи с тем, что при анализе рыночной информации чаще всего приходится иметь дело с нестационарными последовательностями, наибольший интерес может представлять оценка плотности последовательностей малой и средней длины. Это принципиальный момент, который во многом определяет выбор используемого метода оценки.

Для действительно длинных последовательностей, которые содержат больше миллиона значений, с большим успехом могут использоваться гистограммы и P-spline.

Но для последовательности, содержащей 10-20 значений, с эффективным построением гистограмм возникают некоторые проблемы.

Поэтому далее в статье будем ориентироваться в первую очередь на последовательности имеющие длину от 10 до примерно 10000 значений.

1. Методы оценки плотности вероятности

В настоящее время известно большое количество более или менее популярных методов оценки плотности вероятности.

Упоминание о них можно легко найти в сети, воспользовавшись поиском, например, по ключевым словам “оценка плотности вероятности”, “плотность вероятности”,” density estimation” и так далее.

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

Традиционно, для оценки плотности используются гистограммы [1]. Использование гистограмм или сглаженных гистограмм способно обеспечить высокое качество оценок плотности вероятности, но лишь в случае последовательностей большой длины.

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

Поэтому от использования гистограмм пришлось отказаться.

Другим достаточно известным методом оценки является ядерная оценка плотности (kernel density estimation) [2]. Хорошие иллюстрации, поясняющие суть использования ядерного сглаживания можно найти по ссылке [3]. В конечном итоге, несмотря на все недостатки присущие этому методу, выбор пал именно на него. Некоторые моменты, связанные с реализацией этого метода, будут кратко рассмотрены далее.

Среди прочих методов оценки плотности хотелось бы упомянуть о весьма интересном методе, в котором используется алгоритм “Expectation–maximization” [4]. Этот алгоритм позволяет разделить последовательность на отдельные компоненты имеющие, например, нормальное распределение.

После определения параметров отдельных компонент можно просуммировав полученные кривые распределения получить оценку плотности. Упоминание о данном методе можно найти в публикации, размещенной по ссылке [5]. Этот метод, как впрочем, и масса других, при подготовке статьи не был реализован и не тестировался.

Большое количество предлагаемых в различных источниках методов оценки плотности не позволяет практически исследовать каждый из них.

Перейдем к рассмотрению выбранного для реализации метода – метода ядерной оценки плотности.

2. Ядерная оценка плотности вероятности

Ядерная оценка плотности вероятности основывается на методе ядерного сглаживания. Познакомиться с принципами этого метода можно, например, в литературе [6], [7].

Основная идея ядерного сглаживания достаточно проста. Пользователи MetaTrader 5 хорошо знакомы с индикатором Moving Average (MA).

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

Легко заметить, что такое же скользящее окно мы видим и при ядерном сглаживании (например, [3]) только, в этом случае окно симметричное.

Примеры наиболее часто используемых при ядерном сглаживании окон можно найти по ссылке [8]. В случае если при ядерном сглаживании используется регрессия нулевого порядка, то взвешенные значения последовательности, которые попали в окно (ядро), так же как и в случае с MA, просто усредняются.

Такое же применение оконной функции мы можем видеть и в случае, когда сталкиваемся с вопросами связанными с фильтрацией.

Только теперь та же самая процедура будет представлена несколько иначе – при помощи амплитудно-частотной и фазо-частотной характеристики, а ядро (окно) будет называться импульсной характеристикой фильтра.

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

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

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

где

  • x – последовательность длиной n;
  • K – симметричное ядро;
  • h – диапазон, параметр сглаживания.

В дальнейшем для оценок плотности будем использовать только Гауссово ядро:

Как следует из приведенного ранее выражения, плотность в точке X вычисляется как сумма значений ядра для величин, определяемых разностями между значением X и значениями последовательности. При этом точки X, в которых вычисляется плотность, могут и не совпадать со значениями самой последовательности.

Перечислим основные шаги практической реализации алгоритма ядерной оценки плотности.

  1. Производим оценку среднего значения и стандартного отклонения входной последовательности.
  2. Производим нормализацию входной последовательности. Из каждого ее значения вычитаем ранее найденное среднее и делим на величину стандартного отклонения. После такой нормализации исходная последовательность будет иметь нулевое среднее и стандартное отклонение, равное единице. Непосредственно для вычисления плотности такая нормализация не обязательна, но она позволит унифицировать результирующие графики, так как для любой последовательности на шкале X будут располагаться значения, выраженные в единицах стандартного отклонения.
  3. Находим максимальное и минимальное значение в нормализованной последовательности.
  4. Создаем два массива, размер которых должен соответствовать желаемому количеству отображаемых на результирующем графике точек. Например, если график предполагается строить по 200 точкам, то размер массивов должен составлять соответственно по 200 значений.
  5. Оставляем один из созданных массивов для хранения результата, а в другом сформируем значения точек, для которых будет производиться оценка плотности. Для этого в диапазоне между найденными ранее максимальным и минимальным значениями сформируем 200 (для данного примера) равноотстоящих величин и сохраним их в подготовленном массиве.
  6. Используя приведенное ранее выражение, произведем оценку плотности в 200 (для нашего примера) тестовых точках, сохраняя результат в подготовленном на шаге 4 массиве.

Ниже приведена программная реализация этого алгоритма.

#property copyright “2012, victorg” #property link      “https://www.mql5.com” #include class CDens:public CObject   { public:    double            X[];                  int               N;                    double            T[];                  double            Y[];                  int               Np;                  double            Mean;                double            Var;                  double            StDev;                double            H;                 public:    void              CDens(void);    int               Density(double &x[],double hh);    void              NTpoints(int n); private:    void              kdens(double h);   }; void CDens::CDens(void)   {    NTpoints(200);               } void CDens::NTpoints(int n)   {    if(n

Источник: https://www.mql5.com/ru/articles/396

Все о медицине
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: