renuar911's Journal

Глава 23. Мы с Андрюшей постигаем тонкости аппроксимации
Георгий Александров

- Андрюша! Ты сделал домашнее задание на завтра?
- Да, дедуля. Я раскрасил треугольники: равнобедренный, прямоугольный и два простых
- Каких это простых?
- Совсем простых. У них все стороны разные, но у одного все углы острые а другой тупоугольный.
- Ух ты! Много ты знаешь, однако. А я пришел опять по вопросу аппроксимации.
- Деда! Ну сколько можно крутить одну и ту же волынку? Чего мы только с тобой не аппроксимировали!
- Так студенты покоя не дают! Они просят помочь, вот и приходится крутить эту самую... ну, волынку.
- Дедуль! Что-нибудь сложное?
- Наоборот, самое простое, что может быть: линейная регрессия.
- А что такое линейная регрессия?
- Это, внучек, обычная аппроксимация линейной зависимостью.
- Нууу! Нашел чем удивлять студентов! В инете полно мест, где такое делается в режиме онлайн.
- Понимаешь, я обнаружил некоторые тонкости...
- Какие еще тонкости, дедуля? Прямая искривляется?
- Почти попал в точку! Давай не будем ля-ля разводить, а рассмотрим примеры. Вот один шустрик выложил результаты какого-то эксперимента. Первая колонка - это икс, вторая - игрек:

-Деда, я насчитал двадцать четыре точки. Можно я сам попробую наладить твою программу на Yabasic, что методом Монте-Карло?
- Можно, конечно! Я думаю, что ты с легкостью справишься. Ведь столько примеров решали в предыдущих главах.
- Тогда копирую и подстраиваю под нашу задачу:

open #1,"A.txt","r"
open #2,"A1.txt","w"
dim x(100),y(100)
z=.01
for i=1 to 24
input #1 x(i),y(i)
next i
s1=10^150:nn=1000000
a0=1:b0=1
for j=1 to nn
a=a0*(1+z*(ran()-.5))
b=b0*(1+z*(ran()-.5))
s=0
for i=1 to 24
x=x(i):y=y(i)
f=a*x+b
s=s+(y-f)^2
next i
if s<=s1 then
print a,b,s
ak=a:bk=b:sk=s
s1=s
if s<.86 then z=0.00001:fi
a0=a:b0=b
fi
next j
print ak ,bk,sk
print #2,ak using"####.############",bk using"####.############",sk

- Все верно, Андрюша. Теперь запускай и находи коэффициенты линейной регрессии.
- Получилось такое:

y=0.5748∙x+0.7170

Сумма квадратов отклонений равна 0.1236

Теперь можно и строить график. Я это быстренько организую в системе Maple

На рисунке видно, что аппроксимация вполне удовлетворительная. Экспериментальные точки случайно отклоняются от линии регрессии. Что-либо уточнять тут смысла нет. Верно, дедуля?
- У меня нет слов! Ты просто становишься на математически подкованные ноги. Теперь вот второй пример. Его попросил сделать студент на форуме mathhelpplanet.com:

- Деда! Тут абсолютно как первый пример, но количество точек раз, два, ... восемнадцать! То есть программка практически та же, но в двух местах заменить 24 на 18. Ну и файл данных сформировать. Это шестиминутное дело... Все готово! Поехали! Вот и результат линейной аппроксимации:

y = - 0.31094 x + 31.292

Сумма квадратов отклонений равна 1,542. Подкорректирую прогу в Maple и вот он - новый замечательный график:

- Ну что можно, Андрюшенька, сказать? На первый взгляд все удачно. Тенденция ухвачена верно. Но посмотри, - отклонения точек от прямой регрессии не носят случайный характер, а четко обозначает некую синусоиду. Видимо, преподаватель "от фонаря" сочинял координаты так называемых экспериментов и гармонично, циклически искажал потихоньку идеальные точки.
- И что из этого, дедуль? Забракуем задание?
- Нет, зачем? Напротив! Давай представим себе что это не выдумки дядюшки учителя, а самые что ни на есть результаты опытов. И зададимся целью: а можно ли отойти от линейности и улучшить показатель? То есть уменьшить существенно квадратичные отклонения.
- Так это же элементарно! Возьмем полином семнадцатой степени и прямо по точкам им протопаем!
- Только не полином! Мы же когда-то с тобой выяснили, что никакой полином не позволит экстраполировать линию, тем более семнадцатой степени! Да он так рванет куда-нибудь вниз или вверх, что волосы на голове зажужжат.
- Как же нам тогда взволновать функцию?
- А что если, Андрюшенька, нам с тобой аппроксимирующую функцию принять в такой форме:

y = a x + b + c∙sin⁡(k x + m)

- Так тогда придется в прогу включить еще эээ... раз, два, три новых параметра! То есть всего будет пять независимых параметров.
- Да, ну и что? Ты же хотел с полиномом аж восемнадцать параметров вводить! Я в три с лишним раза меньше предлагаю. Короче, ты задание понял и вполне без меня справишься. Мне что-то бабушка твоя сказать хочет. Проблема, значит, появилась на кухне. Сейчас вернусь...
- Деда! Все готово! Вот программа, проверь на всякий пожарный:

open #1,"A.txt","r"
open #2,"A2d.txt","w"
dim x(100),y(100)
z=.001
for i=1 to 18
input #1 x(i),y(i)
next i
s1=10^150:nn=20000000
a0=-.3:b0=30:c0=.3:d0=0.4:k0=0.02
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))
k=k0*(1+z*(ran()-.5))
s=0
for i=1 to 18
x=x(i):y=y(i)
f=a*x+b+c*sin(d*x+k)
s=s+(y-f)^2
next i
if s<=s1 then
print a,b,c,d,k,s
ak=a:bk=b:ck=c:dk=d:kk=k:sk=s
s1=s
if s<1.0 then z=0.00001:fi
a0=a:b0=b:c0=c:d0=d:k0=k
fi
next j
print ak ,bk,ck,dk,kk,sk
print #2,ak using"####.############",bk using"####.############",ck using"####.############",dk using"####.############",kk,sk

- Все вроде четко. Запускай!
- Деда! Смотри - в результате получена значительно более лучшая аппроксимация:

y = -0.31962 x + 31.4487 + 0.35541∙sin⁡( 0.46129 x + 0.0293)

Сумма квадратов отклонений равна всего-навсего 0.541

Сейчас при помощи Maple сформирую рисунок и выясню попутно - насколько улучшилась аппроксимация:

- Вот видишь, внучек, аж на 69 процентов уменьшили сумму отклонений! Это уже серьезное улучшение! Поздравляю тебя с крупным успехом!
- Деда, всего-то 69 процентов? А бывали, случаи, когда точность повышалась в разы?
- Было такое дело. Года три назад ко мне обратился студент. Он привел таблицу каких-то данных... двадцать шесть точек, и попросил помочь произвести линейную аппроксимацию с помощью МНК. Впрочем, сейчас найду эти точки в моей памятке... Так, 25 сентября 2009 года. Вот эти точки:

К сожалению, я не записал результата аппроксимации, но мы с тобой теперь запросто это выполним. Итак, ищем уравнение

y=ax+b

- Деда! Я сам... быстро. Программа такая:

open #1,"A.txt","r"
open #2,"A2.txt","w"
dim x(100),y(100)
z=.01
for i=1 to 26
input #1 x(i),y(i)
next i
s1=10^150:nn=1000000
a0=3:b0=-1
for j=1 to nn
a=a0*(1+z*(ran()-.5))
b=b0*(1+z*(ran()-.5))
s=0
for i=1 to 26
x=x(i):y=y(i)
f=a*x+b
s=s+(y-f)^2
next i
if s<=s1 then
print a,b,s
ak=a:bk=b:sk=s
s1=s
if s<.56 then z=0.00001:fi
a0=a:b0=b
fi
next j
print ak ,bk,sk
print #2,ak,bk,sk

Запускаем ее, минуточку ждем. Вот уравнение! y=2.65097 x -0.665128

Сумма квадратов отклонений 0.558

Теперь строю в Maple сопоставление экспериментов и аппроксимации:

- Теперь, Андрюшенька, анализируем. Вот смотри: вроде все идет гладко, но точки опять же не хаотично разбросаны относительно красной линии. Сначала они резко поднимаются над прямой, образуют бугорок и затем плавно уходят под прямую. Нужно покрупней это исследовать. Давай вот что: проведем прямую между первой и последней точкой. Ты помнишь, как выводить уравнение прямой по двум заданным точкам?
- Да, конечно. Пишу систему двух уравнений, решаю ее и нахожу сначала параметр a=2.722 , а затем и свободный член b=-1.05 . То есть уравнение этой прямой:

y=2.722 x-1.05

График такой:

- Ну а теперь, внучек, даю тебе сложное задание. Слушай внимательно и ничего не перепутай. Твоя задача - усилить получившийся бугорок, которые образуют точки. Нужно сделать следующее: найти отклонение прямой от каждой точки и построить график этих отклонений. Желательно в крупном виде.
- Понял, дедуля. Буду, хоть расколюсь, решать эту задачу... А ты можешь пока свои дела осуществлять. Я тебя позову...
Деда! Вот написал прогу в Maple по результатам вычислений разницы:

restart; with(plots); X := [0, .2, .4, .6, .8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5]; Y := [0, 0.556e-1, .1812, .3168, .4224, .488, .5136, .5092, .4748, .4304, .386, .3316, .2772, .2328, .1886, .154, .1196, 0.952e-1, 0.708e-1, 0.564e-1, 0.42e-1, 0.276e-1, 0.232e-1, 0.88e-2, 0.44e-2, 0]; g1 := plot([`$`([X[i], Y[i]], i = 1 .. 26)], x = 0 .. 5, style = POINT, symbol = BOX, color = black); display(g1);

Запускаю, и - вот он, красавец:

- Отлично, Андрик! Вижу прекрасно, что не ошибся. На что же это похоже? Как ты думаешь?
- Я думаю, что похоже на плотную вероятность.
- Ха-ха! Правильно говорить - на функцию плотности вероятностей. Мы с тобой, помнишь, сколько понастроили? Итак, мой опыт подсказывает: решение нужно искать в виде

- Это я теперь запросто!

open #1,"A.txt","r"
open #2,"A1.txt","w"
dim x(100),y(100)
z=.01
for i=1 to 26
input #1 x(i),y(i)
next i
s1=10^150:nn=1000000
a0=3:b0=-1:c0=3:d0=2:k0=-2
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))
k=k0*(1+z*(ran()-.5))
s=0
for i=1 to 26
x=x(i):y=y(i)
f=a*x+b+c*x^d*exp(k*x)
s=s+(y-f)^2
next i
if s<=s1 then
print a,b,c,d,k,s
ak=a:bk=b:ck=c:dk=d:kk=k:sk=s
s1=s
if s<.0001 then z=0.00001:fi
a0=a:b0=b:c0=c:d0=d:k0=k
fi
next j
print ak ,bk,ck,dk,kk,sk
print #2,ak,bk,ck,dk,kk,sk

Результат прост:

Сумма квадратов отклонений 0.00247 . Вот это да! Ужасно хочется увидеть результат аппроксимации:

restart; with(plots); X := [0, .2, .4, .6, .8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5]; Y := [-1.05, -.45, .22, .9, 1.55, 2.16, 2.73, 3.27, 3.78, 4.28, 4.78, 5.27, 5.76, 6.26, 6.76, 7.27, 7.78, 8.3, 8.82, 9.35, 9.88, 10.41, 10.95, 11.48, 12.02, 12.56]; g1 := plot([`$`([X[i], Y[i]], i = 1 .. 26)], x = 0 .. 5, style = POINT, symbol = BOX, color = black); g2 := plot(2.72*x-1.0542+2.9683*x^2.2526*exp(-1.793*x), x = 0 .. 5, thickness = 2); display(g2, g1);

- Прекрасно, Андрюшенька! Получается, что аппроксимация улучшена почти в 15 раз!
- Ура, дедуля! Вот это я понимаю - улучшение.

Москва
10 января 2012 г.

Current Mood: grateful