Линейная регрессия

[tatkuz1990]

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

[Евгений Машеров]

Во-первых, действительно ли нужна минимизация суммы квадратов отклонений измеренной переменной от оцененной, или же минимизация суммы квадратов отклонений логарифмов (второе на практике встречается куда чаще)?
Во-вторых, есть алгоритмы нелинейной регрессии, основанные на линеаризации, например, Левенберга-Марквардта. В них строится вспомогательная линейная регрессия, в которой регрессоры - производные по параметрам модели, а регрессанд - отклонения наблюдаемой величины Y от расчётной. Коэффициенты регрессии - поправка к значениям коэффициентов нелинейной модели.

[tatkuz1990]

Вот мне бы и хотелось иметь подобную модель с учетом поправок. Мне нужна минимальная сумма квадратов отклонений самой функции, но никак не ее логарифма. Пусть формула будет в три раза большей, но чтобы была явной.
Наверняка в литературе такое имеется.
Про алгоритм Левенберга-Марквардта: http://www.katy.ru/science/LMA.pdf Пришла интересная идея, которую попробую реализовать. Если получится, непременно сообщу. Завтра кровь из носу нужно проблему решить...

[Евгений Машеров]

Формула будет не "в три раза большей", а итеративной. Для модели $y_i=f(x_i,\alpha)$, выбрав начальные значения $\alpha_0$, строим вспомогательную линейную регрессию с регрессандом $u_i=y_i-f(x_i,\alpha^k)$, где k - номер итерации, и регрессорами z_{i,j}=\frac {\partial f(x_j, \alpha^k_j)}{\partial \alpha^k_j}$ вида $u=Z\beta$ и $\beta^k=(Z^TZ)^{-1}Z^Tu$ даёт поправку к альфе.
Повторяем до сходимости. При этом в случае нескольких параметров матрица производных, возможно, будет близка к матрице неполного ранга (мультиколлинеарность), что заставляют применять приёмы повышения обусловленности типа ридж-регрессии (что, собственно, и есть метод Левенберга-Марквардта), для одного параметра такой проблемы нет.

[tatkuz1990]

Метод Левенберга-Марквардта оказался практически малопригодным, поскольку число параметров может доходить до 4-5. Поставленную в теме задачу удалось достаточно просто решить так: в первом приближении рассчитываем два параметра классическим методом линейной регрессии, а затем методом Монте-Карло (с любой заданной точностью) эти два параметра сводим к оптимуму.

[Александрович]

Вот такое нау-хау.
Функция регрессии: $T(t)=T_0 \exp(at)$.
Находим $a$ после логарифмирования.
1. Формируем три параметра $a_{0}=0,9a; a_1=a; a_{2}=1,1a$, то есть 10%-ные отклонения вверх и вниз. Для каждого параметра делаем замену переменных $z_k(t)=\exp(a_kt); k=0,1,2.$ Имеем три функции регрессии вида $T_k(t)=T_{0k}z_k$. По МНК находим $T_{0k}$ и считаем сумму квадратов невязок $S_k(a_k)$. Считая эти 3 пары координатами параболы находим $a'$, при котором парабола имеет минимум.
2. Повторяем п.1 с новым значением $a'$ и с 5%-ным отклонением. Находим $a''$.
3. Повторяем п.1 с новым значением $a''$ и с 1%-ным отклонением. Находим $a'''$.
4. Для функции регрессии $T(t)=T_{0}z'''$ находим $T_0$ по МНК.
Проверил на вашем примере, коэффициенты получаются оптимальные.

[tatkuz1990]

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

[tatkuz1990]

Александрович, я тоже получил оптимальные коэффициенты до 8-го знака после запятой. Какие Ваши результаты?

[Евгений Машеров]

Вообще-то Левенберг-Марквардт как раз для многомерной модели, на 10-20 параметрах работает прекрасно... Да и более справится.

[tatkuz1990]

Монте-Карло работает при неограниченном числе параметров.
Такой оптимум был получен для поверочного примера



Точность до 8-го знака после запятой стабилизировалась.

[Александрович]

tatkuz1990 в сообщении #820838 писал(а):

Александрович, я тоже получил оптимальные коэффициенты до 8-го знака после запятой. Какие Ваши результаты?

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

[Евгений Машеров]

А можно подробнее о Вашей процедуре Монте-Карло?

[ИСН]

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

[tatkuz1990]

Объяснять словами долго, лучше я напишу текст программы и дам комментарии. Итак текст предельно простой:

dim x(100),y(100)
n=4
x(1)=0:x(2)=1:x(3)=2:x(4)=3
y(1)=13.8:y(2)=7.9:y(3)=6.1:y(4)=2.9
rem ЛИНЕЙНАЯ  РЕГРЕССИЯ
for k=1 to n
sx=sx+x(k)
sy=sy+y(k)
sx2=sx2+x(k)^2
slny=slny+log(y(k))
sxlny=sxlny+x(k)*log(y(k))
next k
b=(n*sxlny-sx*slny)/(n*sx2-(sx)^2)
a=slny/n-b/n*sx
rem МЕТОД  МОНТЕ-КАРЛО
z=.00001
a0=a:b0=b
s1=10^150
for j=1 to 3000000
a=a0*(1+z*(ran()-.5))
b=b0*(1+z*(ran()-.5))
s=0
for k=1 to 4
x=x(k):y=y(k)
f=exp(a+b*x)
s=s+(y-f)^2
next k
if s<=s1 then
print a using "##.########",b using "##.########",s using "##.########"
s1=s
a0=a:b0=b:fi
next j

Это элементарный бейсик, тут все понятно, как в блок-схемах.
Сначала составляю массив данных, затем вычисляю коэффициенты линейной регрессии. Далее идет Монте-Карло. Начинаю случайно менять значения двух коэффициентов и для них определять сумму квадратов отклонений. Печатаю новые коэффициенты только тогда, когда предыдущий результат хуже. Число циклов для данной задачи принял 3 млн. Для моего мощного компа - это всего 0,1 с (если не высвечивать на мониторе все вычисления). Если высвечивать (как показано в листинге), то секунд 5 считает. Но это с точностью до 8 знака после запятой. Чем меньше нужна точность, тем естественно быстрее. Но сразу оговорюсь: меня не волнует время решения задачи - пусть хоть неделю пашет. Главное - найти оптимум.
Вот так если коротко. Выше мной приведен самый хвостик расчетов (что на черном фоне).

[ИСН]

Опрокидывать доску - за это бьют канделябрами. Мы-то думали, что Вам нужно решить первоначальную задачу, быстро, многократно.