Сравнительный анализ М-оценок в регрессионных моделях

Таблица 2.7. Результаты оценивания для различных распределений ошибок. Данные матрицы Х распределены по закону Тьюки с параметром зашумления г=0,1 и дисперсиями у12=100, у22=1

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

Таблица 2.8. Результаты дополнительных тестов при ошибках модели, распределенных по закону Тьюки. Данные матрицы Х распределены по закону Тьюки с параметром зашумления г=0,1 и дисперсиями у12=100, у22=1

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

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

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

• Глава 3. Построение модели на реальных данных
• Для построения регрессионной модели на реальных данных был выбран набор данных, включающий уровень пенсии и ВВП на душу населения в 27 странах за 2010 год [2], [14]. Для удобства, таблица с данными вынесена в Приложение № 9.
• В работе строится линейная регрессионная модель зависимости уровня пенсии от ВВП страны на душу населения. Параметры регрессионной модели будут оценены методом наименьших квадратов (МНК) и методом наименьших модулей (МНМ). Результаты работы представлены на Рисунке 3.1 Синими кружками выделены смоделированные данные, желтым цветом выделена прямая, построенная по методу наименьших квадратов, красным - с помощью М-оценок Fair, зеленым - М-оценок Хубера, фиолетовым - МНМ, синим - М-оценок Коши, черным - М-оценки Тьюки, а бирюзовым - М-оценки Вельша.
• Рис. 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. Биржевой лидер [Электронный ресурс].

URL: http://www.profi-forex.org/novosti-dnja/entry1008068487.html - (дата обращения: 11.12.2013).

3. Ботвинкин Е.А., Ранговое оценивание параметров регрессионной модели. [Электронный ресурс].

URL: http://www.hse.ru, - (дата обращения: 29.05.2016).

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

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

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

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

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

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

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

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

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

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

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

14. Статистика стран мира. [Электронный ресурс].

URL: http://iformatsiya.ru - (дата обращения: 11.12.2013).

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

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

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

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

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

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

19. Catherine Stuart. Robust regression. [Электронный ресурс].

URL: maths.dur.ac.uk, - (дата обращения: 17.05.2016).

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

21. John Fox. Robust Regression. Appendix to An R and S-PLUS Companion to Applied Regression. [Электронный ресурс].

URL: www.saedsayad.com/docs/RobustRegression.pdf, - (дата обращения: 27.05.2016).

22. 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

23. Robust Statistics: Theory and Methods. Ricardo A. Maronna, R. Douglas Martin and Vґэctor J. Yohai. 2006 John Wiley & Sons, Ltd ISBN: 0-470-01092-4

24. StatSoft [Электронный ресурс].

URL: http://www.statsoft.ru/home/textbook/glossary/gloss_v.html, - (дата обращения: 17.05.2016).


• Приложение 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.ru
...