Кубическая аппроксимация по методу наименьших квадратов реализуется в Oracle не совсем тривиально.
В предыдущей статье мы обсудили сплайн-аппроксимацию кусочными гиперболами - полиномами 3-го порядка. Мы рассуждали о скользящем окне аппроксимации, но не вдавались в подробности собственно кубической аппроксимации, оценки коэффициентов ak полинома
y = a0 + a1 * x + a2 * x * x + a3 * x * x * x.
Матричную математику по этому вопросу заинтересованный читатель легко найдет по ключевой фразе "полиномиальная аппроксимация методом наименьших квадратов"; мы же сосредоточимся на ее воплощении средствами Oracle Database. И снова, как всегда, заметим, что по нашему мнению база данных является не просто хранилищем данных, а мощным средством их обработки.
Обратимся к аналогии. В Oracle есть пара встроенных агрегатных функций для определения коэффициентов a0, a1 линейной регрессии
y = a0 + a1 * x.
Вообще говоря, это понятно: пользователю довольно часто нужно определить линейный тренд своих фактических зашумленных данных, авиационное - тангаж: идем на взлет или на посадку, нос к небу или к земле, climbing or diving. Эти функции, соответственно,
a0: regr_intercept(y, x)
a1: regr_slope(y, x)
-- говорящие названия: "перехват" и "склон". Обратим внимание, что y и x в них поменяны местами - чтоб потом не запутаться.
Отметим еще раз особенность этих функций - они агрегатные. Это значит, что на вход они принимают таблицу из пар (x, y), а на выход выдают скаляр, число. Это напоминает простые агрегатные функции типа count(), sum(), min(), max(), avg() и т.п.
Мы хотим записать пользовательские агрегатные функции для определения по методу наименьших квадратов коэффициентов кубической регрессии a0, a1, a2, a3. Приступим. Рассмотрим на примере коэффициента a3 - остальные проще.
Функция:
create or replace FUNCTION "REGR_CU_A3" (input number_pair) return number
aggregate using regr_cu_a3_typ
;
Хм... А где же математика? Какая-то декларация, и только. Смотрим по шагам.
Входной параметр имеет какой-то странный тип number_pair - проверим его:
create or replace TYPE "NUMBER_PAIR" as object (x number, y number)
;
Более-менее понятно. А что там стоит после using?
create or replace TYPE "REGR_CU_A3_TYP" as object
(
Sx6 number, Sx5 number, Sx4 number, Sx3 number, Sx2 number, Sx1 number, Sx0 number,
Sx3y number, Sx2y number, Sxy number, Sy number,
static function ODCIAggregateInitialize(sctx in out regr_cu_a3_typ) return number
, member function ODCIAggregateIterate(self in out regr_cu_a3_typ, value in number_pair) return number
, member function ODCIAggregateTerminate(self in regr_cu_a3_typ, returnValue out number, flags in number) return number
, member function ODCIAggregateMerge(self in out regr_cu_a3_typ, ctx2 in regr_cu_a3_typ) return number
)
;
Час от часу не легче. Опять какие-то декларации; что-то объектное, похоже. Как сказал мой сын-школьник по сходному вопросу: "Какая-то инженерная ересь..."
И вот тут-то мы, с полным восторгом, обнаруживаем, что у типа Oracle (что, казалось бы, проще?) может быть - барабанная дробь! - тело:
create or replace TYPE BODY "REGR_CU_A3_TYP" is
static function ODCIAggregateInitialize(sctx IN OUT regr_cu_a3_typ)
return number is
begin
sctx := regr_cu_a3_typ(0,0,0,0,0,0,0,0,0,0,0);
return ODCIConst.Success;
end;
member function ODCIAggregateIterate(self IN OUT regr_cu_a3_typ, value IN number_pair) return number is
begin
-- Смысловые вычисления:
self.Sx6 := self.Sx6 + power(value.x, 6);
self.Sx5 := self.Sx5 + power(value.x, 5);
self.Sx4 := self.Sx4 + power(value.x, 4);
self.Sx3 := self.Sx3 + power(value.x, 3);
self.Sx2 := self.Sx2 + power(value.x, 2);
self.Sx1 := self.Sx1 + value.x;
self.Sx0 := self.Sx0 + 1;
self.Sx3y := self.Sx3y + power(value.x, 3) * value.y;
self.Sx2y := self.Sx2y + power(value.x, 2) * value.y;
self.Sxy := self.Sxy + value.x * value.y;
self.Sy := self.Sy + value.y;
--
return ODCIConst.Success;
end;
member function ODCIAggregateTerminate(self IN regr_cu_a3_typ, returnValue OUT number, flags IN number) return number is
D number := 0;
Da3 number := 0; -- tmp
begin
--возвращаем результат из свойства
D := self.Sx6 * (self.Sx4 * self.Sx2 * self.Sx0 + self.Sx3 * self.Sx1 * self.Sx2 + self.Sx2 * self.Sx3 * self.Sx1 - self.Sx2 * self.Sx2 * self.Sx2 - self.Sx3 * self.Sx3 * self.Sx0 - self.Sx4 * self.Sx1 * self.Sx1)
- self.Sx5 * (self.Sx5 * self.Sx2 * self.Sx0 + self.Sx3 * self.Sx1 * self.Sx3 + self.Sx2 * self.Sx4 * self.Sx1 - self.Sx2 * self.Sx2 * self.Sx3 - self.Sx3 * self.Sx4 * self.Sx0 - self.Sx5 * self.Sx1 * self.Sx1)
+ self.Sx4 * (self.Sx5 * self.Sx3 * self.Sx0 + self.Sx4 * self.Sx1 * self.Sx3 + self.Sx2 * self.Sx4 * self.Sx2 - self.Sx2 * self.Sx3 * self.Sx3 - self.Sx4 * self.Sx4 * self.Sx0 - self.Sx5 * self.Sx1 * self.Sx2)
- self.Sx3 * (self.Sx5 * self.Sx3 * self.Sx1 + self.Sx4 * self.Sx2 * self.Sx3 + self.Sx3 * self.Sx4 * self.Sx2 - self.Sx3 * self.Sx3 * self.Sx3 - self.Sx4 * self.Sx4 * self.Sx1 - self.Sx5 * self.Sx2 * self.Sx2)
;
Da3 := self.Sx3y * (self.Sx4 * self.Sx2 * self.Sx0 + self.Sx3 * self.Sx1 * self.Sx2 + self.Sx2 * self.Sx3 * self.Sx1 - self.Sx2 * self. Sx2 * self. Sx2 - self.Sx3 * self.Sx3 * self.Sx0 - self.Sx4 * self.Sx1 * self.Sx1)
- self.Sx5 * (self.Sx2y * self.Sx2 * self.Sx0 + self.Sx3 * self.Sx1 * self.Sy + self.Sx2 * self.Sxy * self.Sx1 - self.Sx2 * self.Sx2 * self.Sy - self.Sx3 * self.Sxy * self.Sx0 - self.Sx2y * self.Sx1 * self.Sx1)
+ self.Sx4 * (self.Sx2y * self.Sx3 * self.Sx0 + self.Sx4 * self.Sx1 * self.Sy + self.Sx2 * self.Sxy * self.Sx2 - self.Sx2 * self.Sx3 * self.Sy - self.Sx4 * self.Sxy * self.Sx0 - self.Sx2y * self.Sx1 * self.Sx2)
- self.Sx3 * (self.Sx2y * self.Sx3 * self.Sx1 + self.Sx4 * self.Sx2 * self.Sy + self.Sx3 * self.Sxy * self.Sx2 - self.Sx3 * self.Sx3 * self.Sy - self.Sx4 * self.Sxy * self.Sx1 - self.Sx2y * self.Sx2 * self.Sx2)
;
if D != 0 then returnValue := Da3 / D; else returnValue := null; end if;
return ODCIConst.Success;
end;
member function ODCIAggregateMerge(self IN OUT regr_cu_a3_typ, ctx2 IN regr_cu_a3_typ) return number is
begin
--в случае объединения - также перемножаем результаты 2 потоков
--self.mult_value := nvl(self.mult_value,0) * nvl(ctx2.mult_value,0);
return ODCIConst.Success;
end;
end;
Вот тут мы выпадаем в белый творожистый осадок и начинаем думать: что это?
А это и есть решение матричного уравнения (системы уравнений) для определения коэффициента a3 аппроксимирующего полинома 3-го порядка по методу наименьших квадратов в форме пользовательской агрегатной функции. Вычисляем определитель матрицы и ее частный определитель, делим одно на другое. Первый курс технического вуза.
Заметим, кроме a3 есть еще a2, a1, a0 - они попроще, но тоже занимательны.
Отметим, что в данном случае для создания пользовательских агрегатных функций используются встроенные возможности Oracle Database. Я вовсе не уверен, что подобное возможно на других БД.