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

[tatkuz1990]

ewertОчень интересно! Спасибо! И сумма квадратов отклонений 0,00159 , что в три раза меньше, чем у нас...
Сейчас заканчиваю отладку, результаты выложу, возможно даже вечером...

[B@R5uk]

ewert в сообщении #822300 писал(а):

и заодно нормировочный множитель

Так ведь это:

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

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

ewert в сообщении #822300 писал(а):

Если зафиксировать нормировку и варьировать оба параметра, то оптимум достигается на функции $f(x)=x^{8.91207818270}e^{7.8813513381-7.89896643823x}$

Всё именно так с точностью до $10^{-6}$, я вчера -1 у степени забыл. Только у меня сумма квадратов отклонений получилась равна $0.001980059929652$.

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

В рамках данной темы первые же пробы дали такой результат:

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

[ewert]

B@R5uk в сообщении #822324 писал(а):

Всё именно так с точностью до $10^{-6}$,

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

B@R5uk в сообщении #822324 писал(а):

Так ведь это:

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

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

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

Впрочем, вот что получается, если девятку и единичную нормировку фиксировать, а подгонять только один оставшийся параметр:

$f(x)=x^9e^{7.956498273804-7.971180305603x};$

$\sum\big(f(x_i)-y_i\big)^2=2.038603402372\cdot10^{-3};\quad\sigma=1.01\cdot10^{-2};$

Относительные отклонения
   17.54579    0.000016
    0.54254    0.020   
    0.04711    0.23    
    0.01770    0.64    
   -0.01457    1.0     
    0.03249    1.0     
   -0.00055    0.84    
    0.03094    0.55    
    0.00714    0.33    
    0.02471    0.17    
   -0.05215    0.088   
   -0.11752    0.042   
   -0.26341    0.021   
   -0.34896    0.0094  
   -0.48616    0.0045  
   -0.58034    0.0020  
   -0.70588    0.0010  
   -0.79187    0.00048 
   -0.86799    0.00025 
   -0.90333    0.00011 

Как видите -- практически что в лоб, что по лбу.

[tatkuz1990]

B@R5uk в сообщении #822324 писал(а):

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

В рамках данной темы первые же пробы дали такой результат:

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

Вечером дам. Сейчас я не дома, на память ничего не запомнилось - слишком много всего ночью лопатилось... Одно только сейчас могу уверенно говорить: функция распределения принята

$F=1-\exp(ax^{bx...})$

дальше, хоть убейте, не помню... Но еще два параметра задействованы.
Функция плотности распределения получилась, как производная от F. Она довольно длинной оказалась, но для точной математики это роли большой не играет. Главное, что площадь в точности равна единице. Об этом, к моему удивлению, говорил и

Евгений Машеров в сообщении #822247 писал(а):

-- 03 фев 2014, 11:08 --По собственно аппроксимации. Либо исходить из того, что это некая известная функция распределения (гамма, Вейбулла) и подбирать параметры, что в принципе может вывести и на мысли о физическом механизме процесса.
Либо искать аппроксимацию не к функции плотности распределения, а к функции распределения. Что законнее с точки зрения статистики, а с точки зрения техники аппроксимации позволит сразу искать функцию, стремящуюся к нулю при аргументе, стремящемся к минус бесконечности и к единице при аргументе, стремящемся к плюс бесконечности.
Что-то вроде $F(x)=\frac 1 {1+e^{-p(x)}}$ где p(x) монотонно возрастающая функция, стремящуюся к минус бесконечности при аргументе, стремящемся к минус бесконечности и к плюс бесконечности при аргументе, стремящемся к плюс бесконечности.

Но прочитана сия ценная мысль только сегодня днем, когда настал час пробуждения. Видимо, мысли - они материальны и витают, не ведая расстояний.
ewert, насчет "ловли блох" принципиально неверно. Если точно не соблюдается граничное условие, то формула считается лишь интерполяционной. Прогнозы же делать уже рискованно (говорю не об этой конкретной задаче, а в общем случае).

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

Возвращаясь к регрессии.
Различие в результатах модели, сведенной к линейной логарифмированием и с непосредственной минимизацией остатков, связано с тем, что эпсилон в одном случае под экспонентой, а в другом вне её. И логарифмирование действует по-разному. Какой надо применять метод - определяется не статистически, а на уровне содержательного анализа, и чаще полезнее логарифмирование (то, что оно ещё проще - "ну, просто повезло").
Если нужен именно вариант $y=e^{a+bx}+\varepsilon$ и вообще нелинеаризуемая (хотя желательно с достаточно гладкими зависимостями) модель, можно её линеаризовать итеративно.
Находим некое начальное приближение, можно в данном случае взять то, которое получено логарифмированием, и строим вспомогательную линейную модель.
В ней y будут разности между измеренными значениями и вычисленными по текущим приближениям коэффициентов, а x-ами, в смысле регрессорами, будут производные вычисленных значений по параметрам. В данном случае производная по a будет попросту сама экспонента, а по b - она же, умноженная на t. И строим такую вспомогательную регрессию (свободный член не вводим, он лишний), полученные коэффициенты которой будут поправками к коэффициентам исходной нелинейной модели. Потом повторяем и т.п. до сходимости.
2.6319 -0.49
2.614814 -0.475959
2.614640 -0.475879
2.614639 -0.475878
Тут возникает некоторая трудность - регрессоры такой вспомогательной модели будут, как правило, сильно коррелированы, имеет место мультиколлинеарность. Поэтому при расчёте вспомогательной регрессии делается некоторое загрубление (именно его наличием отличается метод Левенберга-Марквардта от описанной "наивной реализации"), но для такой простейшей модели это непринципиально.

[ewert]

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

Если точно не соблюдается граничное условие, то формула считается лишь интерполяционной.

Только не интерполяционной, а кисломолочной; интерполяция тут вовсе не при чём. Граничного условия тут тоже нет.

[tatkuz1990]

Интерполяция в том смысле, что, например, графиком можно пользоваться в интервале от 0 до 4 (в нашем случае). Экстраполировать уже нельзя, так как вдруг кривая рванет в космическое пространство.

[ewert]

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

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

Не рванёт, а экстраполировать и не нужно.

[tatkuz1990]

ewert в сообщении #822381 писал(а):

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

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

Не рванёт, а экстраполировать и не нужно.

В том-то и дело, что бывает нужно. Например, при определении вероятности появления землетрясения в районе строительство АЭС, которое бывает раз в 1000 лет.
Такого времени наблюдений нет, ряд этот во много раз короче и потому функцию обеспеченности экстраполируют. Этим и объясняется необходимость очень тщательной аппроксимации.

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

Обеспеченность: $P=1-\exp \left ( a\, x^{b \, x^c+d} \right )$

Плотность вероятности: $f=\frac{dP}{dx}$

Вечером тщательно рассчитаю коэффициенты ...

[Deggial]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Карантин»
Причина переноса: формулы не оформлены $\TeX$ом

tatkuz1990
Наберите все формулы и термы $\TeX$ом.
Инструкции по оформлению формул здесь или здесь (или в этом видеоролике).
После исправлений сообщите в теме Сообщение в карантине исправлено, и тогда тема будет возвращена.

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»
вернул

[ИСН]

Конструкция $e^{x^x}$ не очень-то физична. А впрочем, четырьмя параметрами можно подогнать что угодно, хоть собаку на Луне.

-- менее минуты назад --

Только у Вас не четыре свободных параметра. Или формула не такая. Смотрите, по физическому смыслу эта P должна быть равна 1 на бесконечности (это при отрицательных $a$ уже выполнено автоматически) и 0 в нуле. А что же это значит? Что $d=0$.

-- менее минуты назад --

А нет, нифига. Это мне показалось, что $d$ входит в первый показатель. А на самом деле - во второй? Тогда всё ОК.

[tatkuz1990]

ИСН, было рассмотрено около десятка конструкций четырехпараметрических формул, но дала результат только приведенная выше. Видимо, не все они беспредельно гибки, а уж собаку даже на земле не станут подгонять.
Теперь все по порядку. Думаю, многим будет полезно. Не знаю, причем тут в первой формуле единица, ибо меня интересует производная, для которой эта единица - что шмель для слона. Мы с коллегой выяснили важное обстоятельство: параметр $a$ должен быть строго отрицательным. Знаки остальных трех параметров пока неизвестны. Чтобы особо не возиться с этим делом, решено было перебрать все варианты знаков. Их оказалось восемь. Была составлена комбинаторная таблица начальных параметров коэффициентов (файл "k8.txt)":

-1 1 1 1
-1 1 1 -1
-1 1 -1 1
-1 -1 1 1
-1 1 -1 -1
-1 -1 1 -1
-1 -1 -1 1
-1 -1 -1 -1 

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

open #1,"c.txt","r"
open #2,"k8.txt","r"
open #3,"8.txt","w"
dim x(100),y(100),f1(100)
z=.01
for i=1 to 20
input #1 x(i),y(i)
next i
for v=1 to 8
print v
input #2 a0,b0,c0,d0
s1=10^100:nn=100000
for j=1 to nn
a=a0*(1+z*(ran()-.5))
b=b0*(1+z*(ran()-.5))
c=c0*(1+z*(ran()-.5))
d=d0*(1+z*(ran()-.5))
s=0
for i=1 to 20
x=x(i)
f1(i)= -a*x^(b*x^c+d-1)*(b*x^c*(c*log(x)+1)+d)*exp(a*x^(b*x^c+d))
s=s+(y(i)-f1(i))^2
next i
if s<=s1 then
ak=a:bk=b:ck=c:dk=d:sk=s
s1=s
a0=a:b0=b:c0=c:d0=d
fi
next j
print #3,ak,bk,ck,dk,sk
next v

Самая длинная строка - это и есть производная, то есть наша аппроксимирующая функция. Методом Монте-Карло по прохождении 100 тыс. циклов печатаются значения $a,b,c,d,s$, где $s - $ сумма квадратов отклонений. Таблица с именем 8.txt получилась такой:

 -0.3740618   0.0616409   0.0120472   3.4532580   0.053564034
 -0.373893   3.6169980   0.0004765  -0.1021912   0.053630939
 -0.324076   3.2485257  -0.5476564   0.9644940   0.001408661
 -0.323134 -10.9611253   0.1314949  15.1276809   0.002662473
 -0.323708   4.7222391  -0.3580727  -0.5164878   0.001503931
-62.496773  -0.1210748   0.1484497  -0.0054732   3.618862915
 -0.394588  -0.0299906  -3.5460164   3.4178397   0.035515651
-54.152601  -0.7514125  -3.7945881  -0.3290997   3.618862915

Третий вариант оказался наиболее близким, так как сумма квадратов отклонений равна всего $0.00141$

Его и рассмотрим более тщательно. Текст проги:

open #1,"c.txt","r"
dim x(100),y(100),f1(100)
z=.00001
for i=1 to 20
input #1 x(i),y(i)
next i
a0=-.32:b0=3.25:c0=-.55:d0=.96
s1=10^100:nn=1000000
for j=1 to nn
a=a0*(1+z*(ran()-.5))
b=b0*(1+z*(ran()-.5))
c=c0*(1+z*(ran()-.5))
d=d0*(1+z*(ran()-.5))
s=0
for i=1 to 20
x=x(i)
f1(i)= -a*x^(b*x^c+d-1)*(b*x^c*(c*log(x)+1)+d)*exp(a*x^(b*x^c+d))
s=s+(y(i)-f1(i))^2
next i
if s<=s1 then
print a,b,c,d,s
s1=s
a0=a:b0=b:c0=c:d0=d
fi
next j

Здесь уже миллион циклов (на моем компе время счета около 5 с.). Скриншот конца итераций получился такой:



Следовательно: $a=-0.323943 \, ; \, b=3.34791\, ; \, c=-0.529493\, ; \, d=0.865155$

Сумма квадратов отклонений $\sum S^2=0.001407$

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

Интеграл аппроксимирующей функции при любых параметрах равен 1 ( см. http://www.wolframalpha.com/input/?i=ev ... ..infty%29 )

[ewert]

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

Сумма квадратов отклонений $\sum S^2=0.001407$

Ну и чем 0.0014 лучше, чем 0.0016 или даже 0.0018? -- правильно: практически ничем.

Если Вы озабочены именно хвостами, то МНК следует в любом случае применять в каком-нибудь весовом варианте. Если не озабочены, то всё это крохоборство с добавлением параметров ровным счётом ничего практически не даёт. Не считая очевидного от него вреда: чем больше параметров вводить, тем сильнее модель будет чувствовать зашумления, которые ей вообще-то положено фильтровать.

[tatkuz1990]

Здесь, возможно, и ничем не отличается, но в других задачах может отличаться. Для математики нет понятия много или мало. Есть min , к которому нужно стремиться, если задача численная. Здесь главное - методическая сторона. Ничего похожего в литературе не повстречалось. По сути - это настоящая нелинейная регрессия безо всяких искажений. Что интересно, можно работать с самыми сложными конструкциями формул с любым количеством параметров. Лично меня и моего шефа это вполне устроило. Значит и кому-нибудь еще пригодится.

ewert! Эта задача больше всего нужна мне и ваша критика мне по барабану.

[tatkuz1990]

.Ой, тут глюк откуда-то. Удалил.