Оценки параметров регрессионных моделей

Рис. 3.1. Зависимость пенсии от уровня ВВП

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



где - вектор с оценками параметров после изменения (символ «о» сверху обозначает выброс - «outlier»), а - вектор с оценками параметров до изменения. Результаты для наглядности представлены в Таблице 3.2, проанализировав которые можно сделать следующий вывод: наилучшим образом на выбросы в данных реагирует М-оценка Тьюки и Вельша, а также М-оценка Хьюбера. Наихудшим - МНК-оценка. Визуализация моделей после изменения данных представлена на Рисунке 3.2

Таблица 3.2 Оценки параметров модели на реальных данных

Вид оценок
МНМ	-41,6041	-12,1811	22,4258	21,6713	866,2821
МНК	-83,3851	-1373,0754	31,2978	106,6865	1668984,53
М-оценки Хьюбера	-68,0623	-79,5903	25,4376	29,8315	152,2011
М-оценки Fair	-69,6524	-143,4568	26,2943	35,3134	5528,4340
М-оценки Коши	-55,0598	-22,6662	23,6606	21,6094	1053,5527
М-оценки Вельша	-19,6959	-14,1779	20,8956	20,6862	30,4922
М-оценки Тьюки	-10,3375	-5,1781	20,3699	20,2362	26,6373

Рис. 3.2 Построение моделей после изменения данных

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

В данной главе был рассмотрен пример применения М-оценок параметров линейной регрессионной модели к реальным данным и была построена линейная модель зависимости уровня пенсии от ВВП страны. Также влияние выбросов в данных на М-, МНК- и МНМ-оценки параметров модели с реальными данными были рассмотрены.

Заключение

В данной работе были рассмотрены М-оценки параметров регрессионной модели, а также изучены их свойства. Были введены различные оценки параметров регрессионной модели, а также реализован алгоритм нахождения М-, МНМ- и МНК-оценок.

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

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

В сравнении с рассматриваемыми конкурентами, МНК наиболее точен для оценивания параметров регрессионной модели с шумами, распределёнными по закону Гаусса, также метод точен при применении к треугольному и «двугорбому» распределению на основе треугольных и гаусcовских величин, однако применение метода с случае, когда распределение ошибок имеет «тяжелые хвосты» не приносит хороших результатов.

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

Таблица 4.1 Рекомендации относительно применения методов оценивания параметров линейной регрессионной модели

Распределение/метод оценивания	МНМ	МНК	М-оценка Хьюбера	М-оценка Fair	М-оценка Коши	М-оценка Вельша	М-оценка Тьюки
Стандартное нормальное распределение
Распределение Лапласа
Распределение Коши
Распределение Стьюдента
Двугорбое распределение как комбинация нормальных
Распределение Тьюки						+	+
Треугольное распределение
Двугорбое распределение как комбинация треугольных
Логистическое распределение

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

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

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

В работе был рассмотрен пример применения М-оценок параметров линейной регрессионной модели к реальным данным и была построена линейная модель зависимости уровня пенсии от ВВП страны. Также влияние выбросов в данных на М-, МНК- и МНМ-оценки параметров модели с реальными данными были рассмотрены. При применении М-оценок к реальным данным, наиболее качественными оценками оказались М-оценки Тьюки и Вельша, также неплохой результат показала М-оценка Хьюбера. Широко распространенный метод наименьших квадратов оказался менее качественным в случае выбросов в реальных данных.

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

Список используемой литературы

1. Белов А. Г., Щедрин Б. М. Относительная эффективность МНК- и МНМ- оценок // "Ломоносовские чтения" (к 300-летию М.В.Ломоносова), 14-23 ноября 2011г. -- М., Макс Пресс Москва, факультет ВМК МГУ, 2011. -- С. 113-114.

2. Вьюгин В.В. Математические основы теории машинного обучения и прогнозирования. - М.: 2013. - 387 с.

3. Гуриев С.М., Колотилин А.Д., Сонин К.И. Цены на нефть и риск национали- зации: о чем говорят панельные данные? // Экономический журнал ВШЭ. 2008. Т. 12, № 2.

4. Дж. Себер «Линейный регрессионный анализ», М.: Мир, 1980, с.54

5. Крыштановский, А. О. Анализ социологических данных с помощью пакета SPSS [Текст]: учеб. пособие для вузов / А. О. Крыштановский; Гос. ун-т -- Высшая школа экономи-ки. -- М. : Изд. дом ГУ ВШЭ, 2006. -- 281, [3] с. -- (Учебники Высшей школы экономики). -- Прил.: с. 225-- 281. -- 2000 экз. -- ISBN 5-7598-0373-5 (в пер.).

6. Линник Ю. В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. -- 2-е изд. -- М., 1962.(математическая теория)

7. Мун С.А. Регрессионный анализ в медико-биологических исследованиях: методические рекомендации/ С.А. Мун, А.Н. Глушков, Т.А. Штернис, С.А. Ларин, С.А. Максимов; ГБОУ ВПО КемГМА Минздравсоцразвития России. - Кемерово: КемГМА, 2012. - 115 с.

8. Научная библиотека [Электронный ресурс].

URL: http://stu.sernam.ru/book_stat1.php?id=1 - (дата обращения: 11.12.2015).

9. Савина, И. А. Методика библиографического описания : практическое пособие . - М. : Либерия-Бибинформ, 2007. - 144 с.

10. Сидоренко Е. В. Методы математической обработки в психологии. -- СПб.: ООО «Речь», 2000. -- 350 с., ил., стр. 240-245.]

11. Смирнов А.В. Анализ финансового состояния коммерческих банков. - М.: Финансы и статистика, 2007. - 225 с.

12. Т. Хеттманспергер "Статистические выводы, основанные на рангах", М.: Финансы и статистика,1987, с. 249

13. Тюрин Ю.Н., Макаров А.А. Анализ данных на компьютере. Изд. 3-е, перераб. и доп./Под ред. В.Э.Фигурнова - М.: ИНФРА-М, 2002. [Электронный ресурс].

URL:http://antonpiter.narod.ru/7361/5semestr/VM_analiz_dannix.PDF - (дата обращения: 29.05.2015).

14. Хьюбер Дж. П. Робастность в статистике: Пер. с англ. - М.: Мир, 1984. 304 с., ил.

15. Bill Jacoby. Regression III: Advanced Methods [Электронный ресурс].

URL:http://polisci.msu.edu/icpsr/regress3, - (дата обращения: 17.05.2016).

16. Encyclopedia of Environments, Vol. 2. Abdel H. El-Shaarawi, Walter W. Piegorsch 2002 John Wiley & Sons, Ltd ISBN: 0 472 89997 6

17. Joiner, Brian L.; Rosenblatt, Joan R. (1971), "Some Properties of the Range in Samples from Tukey's Symmetric Lambda Distributions", Journal of the American Statistical Association 66(334): 394-399

Приложение 1

Визуализация зависимости p-функции Тьюки от b

1) Для реализации алгоритма используется вспомогательная функция

weight_function( ), вычисляющая веса w для различных w-функций М-оценок.

load('random_matrix_X_norm.mat')

X = X(:,1);

c = 4.6851;

b=-10:0.01:10;

error_norm = randn(50,1);

Y = X+error_norm;

Y(1:5) = 100;

X(1:5) = 20;

b_0 = (X'*X)^(-1)*X'*Y;

error = Y - b_0*X;

s = 1.4826*median(abs(error));

for t = 1:size(b,2)

error = Y - b(t)*X;

for k=1:50

x= error(k)/s;

if abs(x) <= c

p(k) = c^2/6*(1-(1-(x/c)^2)^3);

else

p(k) = c^2/6;

end

end

loss_function(t) = sum(p);

end

plot(b,loss_function,'k');

hold on;

error = Y - X;

b_0 = (X'*X)^(-1)*X'*Y; %мнк-оценка первого приближения

eps = 0.0001;

s = 1.4826*median(abs(error));

weight = weight_function(error/s,'tukey');

for i = 1:size(X,1)

W(i,i)=weight(i);

end

b_robust = (X'*W*X)^(-1)*X'*W*Y;

while norm(b_robust-b_0, 2)/norm(b_robust,2) >= eps

error = Y - X*b_robust;

s = 1.4826*median(abs(error));

weight = weight_function(error/s,'tukey');

for i = 1:size(X,1)

W(i,i)=weight(i);

end

b_0 = b_robust;

b_robust = (X'*W*X)^(-1)*X'*W*Y;

end

error = Y - b_robust*X;

s = 1.4826*median(abs(error));

for t=1:50

x= error(t)/s;

if abs(x) <= c

p(t) = c^2/6*(1-(1-(x/c)^2)^3);

else

p(t) = c^2/6;

end

end

loss= sum(p);

plot (b_robust, loss,'*');

xlabel('b');

ylabel ('L');

Приложение 2

Код для вычисления асимптотических относительных эффективностей М-оценок по отношению к МНК- и МНМ-оценкам

Результатом работы скрипта является 6 excel-файлов.

«results_mnk_triug.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНК-оценкам в случае треугольного распределения.

«results_mnm_triug.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНМ-оценкам в случае треугольного распределения.

«results_ass_mnk.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНK-оценкам в случае стандартного нормального распределения, распределения Лапласа, распределения Коши, распределения Стьюдента с 2,3,5,8,13,15 степенями свободы, «двугорбого» распределения на основе двух нормальных, Тьюки с параметрами г=0,1, у12=100, у22=1.

«results_ass_mnm.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНM-оценкам в случае стандартного нормального распределения, распределения Лапласа, распределения Коши, распределения Стьюдента с 2,3,5,8,13,15 степенями свободы, «двугорбого» распределения на основе двух нормальных, Тьюки с параметрами г=0,1, у12=100, у22=1.

«results_mnk_log.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНК-оценкам в случае логистического распределения.

«results_mnm_log.xlsx» содержит АОЭ исследуемых М-оценок по отношению к МНМ-оценкам в случае логистического распределения.

clear all

c_k = 2.3849;

c_h = 1.345;

c_t = 4.6851;

c_w = 2.9846;

c_f = 1.3998;

f_norm = @(u) (1/sqrt(2*pi)*exp(-u.^2/2));

f_lapl = @(u) 1/2*exp(-abs(u));

f_koshi = @(u) 1./(pi*(u.^2+1));

n=2;

f_student_2 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

n=3;

f_student_3 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

n=5;

f_student_5 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

n=8;

f_student_8 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

n=13;

f_student_13 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

n=15;

f_student_15 = @(u) gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1+u.^2/n).^(-(n+1)/2);

f_dv_gaus = @(u) 1/2*1/sqrt(2*pi).*exp(-(u-3).^2./2) + 1/2*1/sqrt(2*pi).*exp(-(u+3).^2./2);

y = 0.1;

s_1 = 100;

s_2 = 1;

f_tukey = @(u) y./(sqrt(s_1*2*pi)).*exp(-(u.^2)./(2*s_1)) + (1-y)./(sqrt(s_2*2*pi)).*exp(-(u.^2)./(2*s_2));

distribution = {@(u) f_norm(u), @(u) f_lapl(u), f_koshi, f_student_2, f_student_3,f_student_5,f_student_8, f_student_13,f_student_15, f_dv_gaus, f_tukey};

v_m_mnk = zeros (5,11);

v_m_mnm = zeros (5,11);

k=1;

%1 - norm, 2 - laplass, 3-koshi

%4 - student with n=2, 5 - student with n=3,

%6 - student with n=5, 7 - student with n=8,

%8 - student with n=13, 9 - student with n=15

%10 - dvug_gauss, 11 - tukey,

for i=1:11

if i==1

s=1;

end

if i==2

s=sqrt(2);

end

if i==3

s=100000; %->inf

end

if i==4

s=10000000; % ->inf

end

if i==5

s=sqrt(3/(3-2));

end

if i==6

s=sqrt(5/(5-2));

end

if i==7

s=sqrt(8/(8-2));

end

if i==8

s=sqrt(13/(13-2));

end

if i==9

s=sqrt(15/(15-2));

end

if i==10

s=sqrt(10);

end

if i==11

s=sqrt(10.9);

end

%HUBER

fun_1 = @(u) (u./s).^2.*distribution{i}(u);

fun_2 = @(u) 2*c_h^2.*distribution{i}(u);

fun_3 = @(u) distribution{i}(u);

v_1 = integral(fun_1,-s*c_h,s*c_h);

v_2 = integral(fun_2,s*c_h,inf);

v_3 = integral(fun_3,-s*c_h,s*c_h);

v_h = (v_1+v_2)/v_3^2;

v_m_mnk(k,i) = 1/v_h;

v_m_mnm(k,i) = 1/(4*s^2*v_h*(distribution{i}(0))^2);

k=k+1;

%FAIR

fun_1 = @(u) ((u./s)./(1+abs((u./s)/c_f))).^2.*distribution{i}(u);

pr_1 = @(u) c_f^2./(abs(u./s) + c_f).^2.*distribution{i}(u);

v_1 = integral(fun_1,-inf,inf);

v_2 = integral(pr_1,-inf,inf);

v_f = (v_1)/v_2^2;

v_m_mnk(k,i) = 1/v_f;

v_m_mnm(k,i) = 1/(4*s^2*v_f*(distribution{i}(0))^2);

k=k+1;

%KOSHI

fun_1 = @(u) ((u./s)./(1+((u./s)/c_k).^2)).^2.*distribution{i}(u);

fun_2 = @(u) ((c_k^2*(c_k^2-(u./s).^2)./(c_k^2+(u./s).^2).^2)).*distribution{i}(u);

v_k = integral(fun_1,-inf,inf)/integral(fun_2,-inf,inf)^2;

v_m_mnk(k,i) = 1/v_k;

v_m_mnm(k,i) = 1/(4*s^2*v_k*(distribution{i}(0))^2);

k=k+1;

%WELSH

fun_1 = @(u) ((u./s).*exp(-((u./s)./c_w).^2)).^2.*distribution{i}(u);

pr_1 = @(u) exp(-(u./s).^2./c_w^2).*(c_w^2-2*(u./s).^2)./c_w.^2.*distribution{i}(u);

v_1 = integral(fun_1,-inf,inf);

v_2 = integral(pr_1,-inf,inf);

v_w = (v_1)/v_2^2;

v_m_mnk(k,i) = 1/v_w;

v_m_mnm(k,i) = 1/(4*s^2*v_w*(distribution{i}(0))^2);

k=k+1;

%TUKEY

fun_1 = @(u) ((u./s).*(1-(u./(s*c_t)).^2).^2).^2.*distribution{i}(u);;

pr_1 = @(u) (5*(u./s).^4/c_t^4 - 6*(u./s).^2/c_t^2 +1).*distribution{i}(u);;

v_1 = integral(fun_1,-s*c_t,s*c_t);

v_2 = integral(pr_1,-s*c_t,s*c_t);

v_t = (v_1)/v_2^2;

v_m_mnk(k,i) = 1/v_t;

v_m_mnm(k,i) = 1/(4*s^2*v_t*(distribution{i}(0))^2);

k=1;

end

v_m_mnk = v_m_mnk';

v_m_mnm = v_m_mnm';

filename = 'results_ass_mnk.xlsx';

xlswrite(filename,v_m_mnk,1,'A1:E11')

filename = 'results_ass_mnm.xlsx';

xlswrite(filename,v_m_mnm,1,'A1:E11')

%logistick

f_logist = @(u) exp(-u)./(1+exp(-u)).^2;

s = sqrt(pi^2/3);

%HUBER

fun_1 = @(u) (u./s).^2.*f_logist(u);

fun_2 = @(u) 2*c_h^2.*f_logist(u);

fun_3 = @(u) f_logist(u);

v_1 = integral(fun_1,-s*c_h,s*c_h);

v_2 = integral(fun_2,s*c_h,100);

v_3 = integral(fun_3,-s*c_h,s*c_h);

v_h = (v_1+v_2)/v_3^2;

mnk_log(1) = 1/v_h;

mnm_log(1) = 1/(4*s^2*v_h*(f_logist(0))^2);

%FAIR

fun_1 = @(u) ((u./s)./(1+abs((u./s)/c_f))).^2.*f_logist(u);

pr_1 = @(u) c_f^2./(abs(u./s) + c_f).^2.*f_logist(u);

v_1 = integral(fun_1,-100,100);

v_2 = integral(pr_1,-100,100);

v_f = (v_1)/v_2^2;

mnk_log(2) = 1/v_f;

mnm_log(2) = 1/(4*s^2*v_f*(f_logist(0))^2);

%KOSHI

fun_1 = @(u) ((u./s)./(1+((u./s)/c_k).^2)).^2.*f_logist(u);

fun_2 = @(u) ((c_k^2*(c_k^2-(u./s).^2)./(c_k^2+(u./s).^2).^2)).*f_logist(u);

v_k = integral(fun_1,-100,100)/integral(fun_2,-100,100)^2;

mnk_log(3) = 1/v_k;

mnm_log(3) = 1/(4*s^2*v_k*(f_logist(0))^2);

%WELSH

fun_1 = @(u) ((u./s).*exp(-((u./s)./c_w).^2)).^2.*f_logist(u);

pr_1 = @(u) exp(-(u./s).^2./c_w^2).*(c_w^2-2*(u./s).^2)./c_w.^2.*f_logist(u);

v_1 = integral(fun_1,-100,100);

v_2 = integral(pr_1,-100,100);

v_w = (v_1)/v_2^2;

mnk_log(4) = 1/v_w;

mnm_log(4) = 1/(4*s^2*v_w*(f_logist(0))^2);

%TUKEY

fun_1 = @(u) ((u./s).*(1-(u./(s*c_t)).^2).^2).^2.*f_logist(u);

pr_1 = @(u) (5*(u./s).^4/c_t^4 - 6*(u./s).^2/c_t^2 +1).*f_logist(u);

v_1 = integral(fun_1,-s*c_t,s*c_t);

v_2 = integral(pr_1,-s*c_t,s*c_t);

v_t = (v_1)/v_2^2;

mnk_log(5) = 1/v_t;

mnm_log(5) = 1/(4*s^2*v_t*(f_logist(0))^2);

filename = 'results_mnk_log.xlsx';

xlswrite(filename,mnk_log,1,'A1:E1')

filename = 'results_mnm_log.xlsx';

xlswrite(filename,mnm_log,1,'A1:E1')

%TRIUGOLNOE

a=-2;

b=2;

f_triug = @(u) 2/(b-a) - 2/(b-a)^2.*abs(a+b-2.*u);

s=sqrt((b-a)^2/24);

%HUBER

fun_1 = @(u) (u./s).^2.*f_triug(u);

fun_2 = @(u) 2*c_h^2.*f_triug(u);

fun_3 = @(u) f_triug(u);

v_1 = integral(fun_1,-s*c_h,s*c_h);

v_2 = integral(fun_2,s*c_h,b);

v_3 = integral(fun_3,-s*c_h,s*c_h);

v_h = (v_1+v_2)/v_3^2;

mnk_triug(1) = 1/v_h;

mnm_triug(1) = 1/(4*s^2*v_h*(f_triug(0))^2);

%FAIR

fun_1 = @(u) ((u./s)./(1+abs((u./s)/c_f))).^2.*f_triug(u);

pr_1 = @(u) c_f^2./(abs(u./s) + c_f).^2.*f_triug(u);

v_1 = integral(fun_1,a,b);

v_2 = integral(pr_1,a,b);

v_f = (v_1)/v_2^2;

mnk_triug(2) = 1/v_f;

mnm_triug(2) = 1/(4*s^2*v_f*(f_triug(0))^2);

%KOSHI

fun_1 = @(u) ((u./s)./(1+((u./s)/c_k).^2)).^2.*f_triug(u);

fun_2 = @(u) ((c_k^2*(c_k^2-(u./s).^2)./(c_k^2+(u./s).^2).^2)).*f_triug(u);

v_k = integral(fun_1,a,b)/integral(fun_2,a,b)^2;

mnk_triug(3) = 1/v_k;

mnm_triug(3) = 1/(4*s^2*v_k*(f_triug(0))^2);

%WELSH

fun_1 = @(u) ((u./s).*exp(-((u./s)./c_w).^2)).^2.*f_triug(u);

pr_1 = @(u) exp(-(u./s).^2./c_w^2).*(c_w^2-2*(u./s).^2)./c_w.^2.*f_triug(u);

v_1 = integral(fun_1,a,b);

v_2 = integral(pr_1,a,b);

v_w = (v_1)/v_2^2;

mnk_triug(4) = 1/v_w;

mnm_triug(4) = 1/(4*s^2*v_w*(f_triug(0))^2);

%TUKEY

fun_1 = @(u) ((u./s).*(1-(u./(s*c_t)).^2).^2).^2.*f_triug(u);

pr_1 = @(u) (5*(u./s).^4/c_t^4 - 6*(u./s).^2/c_t^2 +1).*f_triug(u);

v_1 = integral(fun_1,a,b);

v_2 = integral(pr_1,a,b);

v_t = (v_1)/v_2^2;

mnk_triug(5) = 1/v_t;

mnm_triug(5) = 1/(4*s^2*v_t*(f_triug(0))^2);

filename = 'results_mnk_triug.xlsx';

xlswrite(filename,mnk_triug,1,'A1:E1')

filename = 'results_mnm_triug.xlsx';

xlswrite(filename,mnm_triug,1,'A1:E1')

Приложение 3

Код для визуализации функций p, ф, w

1) Используется дополнительная функция plot_3_graph( ), позволяющая вывести в одном окне 3 графика.

function [ smt ] = plot_3_graph(x, ro, psi, w )

smt = 1;

subplot(1,3,1);

plot (x, ro);

ylabel('\rho(u)');

xlabel ('u');

subplot(1,3,2);

plot (x, psi);

ylabel ('\psi(u)');

xlabel ('u');

subplot(1,3,3);

plot (x, w);

ylabel ('w(u)');

xlabel ('u');

end

2) Визуализации функций

x = -10:0.0001:10;

ro_huber = zeros (size(x));

phi_huber = zeros (size(x));

w_huber = zeros (size(x));

c = 1.345;

for k = 1:size(x,2)

if abs(x(k)) <= c

ro_huber(k) = (x(k)^2)/2;

phi_huber(k) = x(k);

w_huber(k) = 1;

else

ro_huber(k) = c*(abs(x(k)) - c/2);

phi_huber(k)=c*sign(x(k));

w_huber(k) = c/abs(x(k));

end

end

plot_3_graph (x, ro_huber, phi_huber, w_huber);

figure;

ro_fair = zeros (size(x));

phi_fair = zeros (size(x));

w_fair = zeros (size(x));

c = 1.3998;

for k = 1:size(x,2)

ro_fair(k) = c^2*(abs(x(k))/c - log(1+abs(x(k)/c)));

phi_fair(k) = x(k)/(1+abs(x(k))/c);

w_fair(k) = 1/(1+abs(x(k))/c);

end

plot_3_graph (x, ro_fair, phi_fair, w_fair);

figure;

ro_koshi = zeros (size(x));

phi_koshi= zeros (size(x));

w_koshi = zeros (size(x));

c = 2.3849;

for k = 1:size(x,2)

ro_koshi(k) = c^2/2*log(1+(x(k)/c)^2);

phi_koshi(k) = x(k)/(1+((x(k))/c)^2);

w_koshi(k) = 1/(1+((x(k))/c)^2);

end

plot_3_graph (x, ro_koshi, phi_koshi, w_koshi);

figure;

ro_g = zeros (size(x));

phi_g= zeros (size(x));

w_g = zeros (size(x));

for k = 1:size(x,2)

ro_g(k) = x(k)^2/2/(1+x(k)^2);

phi_g(k) = x(k)/(1+(x(k))^2)^2;

w_g(k) = 1/(1+(x(k))^2)^2;

end

plot_3_graph (x, ro_g, phi_g, w_g);

figure;

ro_welsh = zeros (size(x));

phi_welsh= zeros (size(x));

w_welsh = zeros (size(x));

c = 2.9846;

for k = 1:size(x,2)

ro_welsh(k) = c^2/2*(1-exp(-(x(k)/c)^2));

phi_welsh(k) = x(k) * exp(-(x(k)/c)^2);

w_welsh(k) = exp(-(x(k)/c)^2);

end

plot_3_graph (x, ro_welsh, phi_welsh, w_welsh);

figure;

ro_tukey = zeros (size(x));

phi_tukey = zeros (size(x));

w_tukey = zeros (size(x));

c = 4.6851;

for k = 1:size(x,2)

if abs(x(k)) <= c

ro_tukey(k) = c^2/6*(1-(1-(x(k)/c)^2)^3);

phi_tukey(k) = x(k)*(1-(x(k)/c)^2)^2;

w_tukey(k) = (1-(x(k)/c)^2)^2;

else

ro_tukey(k) = c^2/6;

phi_tukey(k)=0;

w_tukey(k) = 0;

end

end

plot_3_graph (x, ro_tukey, phi_tukey, w_tukey);

Приложение 4

Пример работы М-оценок

p = randn(1,10);

x = 0:1:9;

y = x + p;

x = [x 15];

y = [y 1];

scatter(x,y, 'filled');grid on;

hold on;

types = cellstr(['lad ';'ols '; 'huber';'fair ';'koshi'; 'welsh';'tukey']);

types_color = (['m'; 'y';'g';'r';'b'; 'c';'k']);

x = [ones(size(x,2),1) x'];

for k = 1:size(types,1)

b = IRLS(x,y,'ols',types(k));

plot (x, b(2)*x + b(1), types_color(k));

hold on;

err(:,k) = y'-x*b;

end

xlabel ('x');

grid on;

ylabel ('y');

Приложение 5

Алгоритм для вычисления IRLS

1) Для реализации алгоритма используется вспомогательная функция weight_function( ), вычисляющая веса w для различных w-функций М-оценок.

function [ weight ] = weight_function(x,f_name)

if strcmp(f_name, 'huber') == 1

c = 1.345;

weight = zeros(size(x));

for k=1:size (weight,1)

if abs(x(k)) <= c

weight(k) = 1;

else

weight(k) = c/abs(x(k));

end

end

end

if strcmp(f_name, 'fair') == 1

c = 1.3998;

weight = zeros(size(x));

for k=1:size (weight,1)

weight(k) = 1/(1+abs(x(k))/c);

end

end

if strcmp(f_name, 'koshi') == 1

c = 2.3849;

weight = zeros(size(x));

for k=1:size (weight,1)

weight(k) = 1/(1+((x(k))/c)^2);

end

end

if strcmp(f_name, 'welsh') == 1

c = 2.9846;

weight = zeros(size(x));

for k=1:size (weight,1)

weight(k)=exp(-(x(k)/c)^2);

end

end

if strcmp(f_name, 'tukey') == 1

c = 4.6851;

weight = zeros(size(x));

for k=1:size (weight,1)

if abs(x(k)) <= c

weight(k) = (1-(x(k)/c)^2)^2;

else

weight(k) = 0;

end

end

end

end

2) Функция для вычисления IRLS.

function [ b_robust ] = IRLS(x,y, initial_b_type,weighted_f_name)

X = x;

if strcmp(weighted_f_name,'ols') == 1 || strcmp(weighted_f_name,'lad') == 1

if strcmp(weighted_f_name,'ols') == 1

b_robust = (X'*X)^(-1)*X'*y';

else

[b_robust, fval]=fminsearch(@(b_robust)sum(abs(y'-X*b_robust)), ones(size(X,2),1));

end

else

if strcmp(initial_b_type,'ols') == 1

b_0 = (X'*X)^(-1)*X'*y';

end

eps = 0.0001;

error = y' - X*b_0;

s = 1.4826*median(abs(error));

weight = weight_function(error/s,weighted_f_name);

for i = 1:size(X,1)

W(i,i)=weight(i);

end

b_robust = (X'*W*X)^(-1)*X'*W*y';

while norm(b_robust-b_0, 2)/norm(b_robust,2) >= eps

error = y' - X*b_robust;

s = 1.4826*median(abs(error));

weight = weight_function(error/s,weighted_f_name);

for i = 1:size(X,1)

W(i,i)=weight(i);

end

b_0 = b_robust;

b_robust = (X'*W*X)^(-1)*X'*W*y';

end

end

end

Приложение 6

Визуализация функций плотностей распределения Гаусса и Тьюки

x = -8:0.001:8;

Y = normpdf(x,0,1);%norm

plot(x,Y)

xlabel('x');

ylabel('f(x)');

hold on;

Y = 0.1*normpdf(x,0,10) + (0.9)*normpdf(x,0,1); %tukey

plot(x,Y,'r')

xlabel('x');

ylabel('f(x)');

Приложение 7

Результаты вычисления М-оценок, МНК- и МНМ-оценок

1) Для реализации используется функция IRLS() [см. Приложение 5]

2) Результаты вычисления записываются в отдельный excel-файл «results.xlsx», содержащий 12 строк и 7 столбцов, соответствующий различным распределениям и оценкам.

clear all;

n=50;

m=2;

n_distributions = 12;

b = [3;-7;8];

b_robust_vector = zeros(m+1,1);

X = 5*rand(n,m);

X = [ones(n,1) X];

types = cellstr(['lad ';'ols '; 'huber';'fair ';'koshi'; 'welsh';'tukey']);

kriteriy = zeros(size(types,1),n_distributions);

b_robust_mean = zeros(m+1, size(types,1));

error_matrix = zeros(n,n_distributions);rng(1);

for dist_column = 1:n_distributions

for t = 1:size(types,1)

for k=1:2000

%1norn

error_norm = randn(n,1);

%2laplass

z1=-log(rand(n,1));

z2=-log(rand(n,1));

error_laplass=z1-z2;

%3koshi

error_koshi=tan(pi*rand(n,1) - pi/2);

%4,5,6,7student

error_stud_2 = tinv(rand(n,1),2);

error_stud_3 = tinv(rand(n,1),3);

error_stud_8 = tinv(rand(n,1),8);

error_stud_15 = tinv(rand(n,1),15);

%8dvugorbovoe

p=rand(n,1);

error_gaus_dvud=randn(n,1)+3;

eps=randn(n,1)-3;

I=find(p>=0.5);

error_gaus_dvud(I)=eps(I);

%9tukey

p=rand(n,1);

error_tukey=randn(n,1);

eps=sqrt(100)*randn(n,1);

I=find(p>=0.9);

error_tukey(I)=eps(I);

%10triangle

u1=2*rand(n,1)-1;

u2=2*rand(n,1)-1;

error_triangle=u1+u2;

%11double_triangle

p=rand(n,1);

u1=rand(n,1);

u2=rand(n,1);

error_doub_triang=u1+u2;

v1=-rand(n,1);

v2=-rand(n,1);

eps=v1+v2;

I=find(p>=0.5);

error_doub_triang(I)=eps(I);

%12logistick

error_logistick=-log((1./rand(n,1)) -1);

error_matrix = [error_norm error_laplass error_koshi error_stud_2 error_stud_3 error_stud_8 error_stud_15];

error_matrix = [error_matrix error_gaus_dvud error_tukey error_triangle error_doub_triang error_logistick];

Y = X*b+error_matrix(:,dist_column);

b_robust = IRLS(X,Y','ols',types(t));

b_robust_vector(:,k) = b_robust;

end

for p = 1:size(b_robust_vector,2)

d_vector_error(1,p) = sum((b_robust_vector(:,p) - b).^2);

end

avg_error_b(t,dist_column) = mean(d_vector_error');

end

end

avg_error_b = avg_error_b';

filename = 'results.xlsx';

xlswrite(filename,avg_error_b,1,'A1:G12')

Приложение 8

Дополнительные тесты для распределения Тьюки

1) Для реализации используется функция IRLS() [см. Приложение 5]

2) Дополнительны тесты для Распределения Тьюки. Матрица X может содержать данные, распределенные нормально ('random_matrix_X_norm), по закону Коши ('random_matrix_X_koshi'), либо по закону Тьюки ('random_matrix_X_tukey'). Результат программы - excel-файл 'results_tukey.xlsx';

clear all

load('random_matrix_X_norm.mat')

n=50;

m=2;

types = cellstr(['lad ';'ols '; 'huber';'fair ';'koshi'; 'welsh';'tukey']);

X = [ones(n,1) X];

b = [3;-7;8];

b_robust_vector = zeros(m+1,1);

param = [0.05 0.1 0.15 0.2 0.3];

for r=1:9

for t = 1:size(types,1)

for k=1:2000

p=rand(n,1);

disp_1 = disp1(r);

disp_2=disp2(r);

error_tukey=sqrt(1)*randn(n,1);

eps=sqrt(100)*randn(n,1);

I=find(p>=1-param(r));

error_tukey(I)=eps(I);

Y = X*b+error_tukey;

b_robust = IRLS(X,Y','ols',types(t));

b_robust_vector(:,k) = b_robust;

end

for p = 1:size(b_robust_vector,2)

d_vector_error(1,p) = sum((b_robust_vector(:,p) - b).^2);

end

avg_error_b(t,r) = mean(d_vector_error');

end

end

avg_error_b = avg_error_b';

filename = 'results_tukey.xlsx';

xlswrite(filename,avg_error_b,1,'A1:G9')

Приложение 9

Данные по 27 странам. Уровень пенсии и ВВП

№	Страна	Пенсия, USD	ВВП, млн. USD
1	Дания	2800	37,00
2	Финляндия	1900	35,30
3	Норвегия	1542	59,10
4	Израиль	1350	29,50
5	Германия	1200	35,90
6	Испания	1190	29,50
7	США	1164	47,40
8	Швейцария	874	42,90
9	Швеция	833	39,00
10	Япония	717	34,20
11	Великобритания	700	36,12
12	Франция	700	33,30
13	Канада	667	39,60
14	Италия	583	30,70
15	Венгрия	400	19,00
16	Польша	380	18,80
17	Литва	298	15,90
18	Россия	285	15,90
19	Болгария	280	12,80
20	Казахстан	210	12,50
21	Азербайджан	202	11,00
22	Белоруссия	175	7,41
23	Украина	142	6,70
24	Аргентина	96	14,70
25	Молдова	80	2,50
26	Узбекистан	55	3,10
27	Грузия	40	4,80

Приложение 10

Построение модели зависимости пенсии от уровня ВВП

1) Для реализации используется функция IRLS() [см. Приложение 5]

2) Построение модели

load('pension_and_gdp.mat')

p=fix(rand(1)*5);

Y = X(:,1);

Y(p) = Y(p)*5;

p=fix(rand(1)*5);

Y(p) = Y(p)*5;

X=[ones(27,1) X(:,2)];

scatter(X(:,2),Y, 'filled');grid on;

hold on;

types = cellstr(['lad ';'ols '; 'huber';'fair ';'koshi'; 'welsh';'tukey']);

types_color = (['m'; 'y';'g';'r';'b'; 'c';'k']);

for t=1:size(types,1)

b_robust = IRLS(X,Y','ols',types(t));

b_robust_vector(:,t) = b_robust;

plot(X(:,2), b_robust(2)*X(:,2) + b_robust(1), types_color(t));

hold on;

end

xlabel('ВВП, млн долл. США');

ylabel('Пенсия, долл. США');

Размещено на Allbest.ur

...