Сравните устойчивые методы регрессии

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

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

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

Моделируйте данные

Сгенерируйте случайную выборку 100 ответов из линейной модели

$y_{t} = 1 + 2 x_{t} + ε_{t}$

где:

$x$ вектор равномерно расположенных с интервалами значений от 0 до 2.
$ε_{t} \sim N (0, 0 . 5^{2}) .$

rng('default');
n = 100;
x = linspace(0,2,n)';
b0 = 1;
b1 = 2;
sigma = 0.5;
e = randn(n,1);
y = b0 + b1*x + sigma*e;

Создайте влиятельные выбросы путем расширения всего соответствия ответов $x < 0 . 25$ фактором 3.

y(x < 0.25) = y(x < 0.25)*3;

Отобразите данные на графике. Сохраните график для дальнейших графиков.

figure;
plot(x,y,'o');
h = gca;
xlim = h.XLim';
hl = legend('Data');
xlabel('x');
ylabel('y');
title('Regression Techniques Comparison')
hold on;

Оцените линейную модель

Оцените коэффициенты и ошибочное отклонение при помощи простой линейной регрессии. Постройте линию регрессии.

LSMdl = fitlm(x,y)

LSMdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate      SE       tStat       pValue  
                   ________    _______    ______    __________

    (Intercept)     2.6814     0.28433    9.4304    2.0859e-15
    x1             0.78974     0.24562    3.2153     0.0017653


Number of observations: 100, Error degrees of freedom: 98
Root Mean Squared Error: 1.43
R-squared: 0.0954,  Adjusted R-Squared: 0.0862
F-statistic vs. constant model: 10.3, p-value = 0.00177

plot(xlim,[ones(2,1) xlim]*LSMdl.Coefficients.Estimate,'LineWidth',2);
hl.String{2} = 'Least Squares';

LSMdl является подходящим объектом модели LinearModel. Прерывание и наклон, кажется, соответственно выше и ниже, чем они должны быть. Линия регрессии может сделать плохие прогнозы для любого $x < 1$ и $x > 1.6$ .

Оцените байесовую модель линейной регрессии с рассеянным предшествующим распределением

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

PriorDiffuseMdl = bayeslm(1);

PriorDiffuseMdl является объектом модели diffuseblm, который характеризует объединенное предшествующее распределение.

Оцените следующую из Байесовой модели линейной регрессии. Постройте линию регрессии.

PosteriorDiffuseMdl = estimate(PriorDiffuseMdl,x,y);

Method: Analytic posterior distributions
Number of observations: 100
Number of predictors:   2
 
           |  Mean     Std         CI95        Positive      Distribution     
------------------------------------------------------------------------------
 Intercept | 2.6814  0.2873  [ 2.117,  3.246]    1.000   t (2.68, 0.28^2, 98) 
 Beta      | 0.7897  0.2482  [ 0.302,  1.277]    0.999   t (0.79, 0.25^2, 98) 
 Sigma2    | 2.0943  0.3055  [ 1.580,  2.773]    1.000   IG(49.00, 0.0099)

plot(xlim,[ones(2,1) xlim]*PosteriorDiffuseMdl.Mu,'--','LineWidth',2);
hl.String{3} = 'Bayesian Diffuse';

PosteriorDiffuseMdl является объектом модели conjugateblm, который характеризует объединенное апостериорное распределение линейных параметров модели. Оценки Байесовой модели линейной регрессии с предшествующим рассеянным почти равны тем из простой модели линейной регрессии. Обе модели представляют наивный подход к влиятельным выбросам, то есть, методы обрабатывают выбросы как любое другое наблюдение.

Оцените модель регрессии с ошибками ARIMA

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

regMdl = regARIMA(0,0,0);
regMdl.Distribution = struct('Name','t','DoF',3);

regMdl является объектом модели regARIMA. Это - шаблон для оценки.

Оцените модель регрессии с ошибками ARIMA. Постройте линию регрессии.

estRegMdl = estimate(regMdl,y,'X',x);

 
    Regression with ARMA(0,0) Error Model (t Distribution):
 
                  Value     StandardError    TStatistic      PValue  
                 _______    _____________    __________    __________

    Intercept     1.4613         0.154         9.4892       2.328e-21
    Beta(1)       1.6531       0.12939         12.776      2.2245e-37
    Variance     0.93116        0.1716         5.4263      5.7546e-08
    DoF                3             0            Inf               0

plot(xlim,[ones(2,1) xlim]*[estRegMdl.Intercept; estRegMdl.Beta],...
    'LineWidth',2);
hl.String{4} = 'regARIMA';

estRegMdl является объектом модели regARIMA, содержащим результаты оценки. Поскольку t распределение - больше рассеянное, линия регрессии приписывает больше изменчивости влиятельным выбросам, чем к другим наблюдениям. Поэтому линия регрессии, кажется, лучшая прогнозирующая модель, чем другие модели.

Реализуйте регрессию квантиля Используя мешок деревьев регрессии

Вырастите мешок 100 деревьев регрессии. Задайте 20 для минимального листового размера.

QRMdl = TreeBagger(100,x,y,'Method','regression','MinLeafSize',20);

QRMdl является подходящим объектом модели TreeBagger.

Предскажите средние ответы для всех наблюдаемых $x$ значения, то есть, регрессия квантиля реализации. Постройте прогнозы.

qrPred = quantilePredict(QRMdl,x);
plot(x,qrPred,'LineWidth',2);
hl.String{5} = 'Quantile';
hold off;

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

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

Реализуйте устойчивую байесовую линейную регрессию

Рассмотрите Байесовую модель линейной регрессии, содержащую один предиктор, t распределил отклонение воздействия с профилируемым параметром степеней свободы $ν$ Пусть:

$λ_{j} \sim I G (ν / 2, 2 / ν)$ .
$ε_{j} | λ_{j} \sim N (0, λ_{j} σ^{2})$
$f (β, σ^{2}) \propto \frac{1}{σ^{2}}$

Эти предположения подразумевают:

$ε_{j} \sim t (0, σ^{2}, ν)$
$λ_{j} | ε_{j} \sim I G (\frac{ν + 1}{2}, \frac{2}{ν + ε_{j}^{2} / σ^{2}})$

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

Задайте сетку значений для $ν$ .

1 соответствует распределению Коши.
2.1 означает, что среднее значение четко определено.
4,1 средних значения, что отклонение четко определено.
100 означает, что распределение приблизительно нормально.

nu = [0.01 0.1 1 2.1 5 10 100];
numNu = numel(nu);

Для этой проблемы сэмплер Гиббса является подходящим, чтобы оценить коэффициенты, потому что можно моделировать параметры Байесовой модели линейной регрессии, обусловленной на $λ$ , и затем моделируйте $λ$ от его условного распределения.

Реализуйте этот сэмплер Гиббса.

Чертите параметры от апостериорного распределения $β, σ^{2} | y, x, λ$ . Выкачайте наблюдения $λ$ , создайте рассеянную предшествующую модель с двумя коэффициентами регрессии и чертите набор параметров от следующего. Первый коэффициент регрессии соответствует прерыванию, поэтому укажите, что bayeslm не включает тот.
Вычислите невязки.
Чертите значения от условного выражения, следующего из $λ$ .

Для каждого значения $ν$ , запустите сэмплер Гиббса для 20 000 итераций и примените электротермотренировку 5 000. Предварительно выделите для следующих ничьих и инициализируйте $λ$ к вектору из единиц.

rng(1);
m = 20000;
burnin = 5000;
lambda = ones(n,m + 1,numNu); % Preallocation
estBeta = zeros(2,m + 1,numNu);
estSigma2 = zeros(1,m + 1,numNu);

% Create diffuse prior model.
PriorMdl = bayeslm(2,'Model','diffuse','Intercept',false);

for p = 1:numNu
    for j = 1:m
        % Scale observations.
        yDef = y./sqrt(lambda(:,j,p)); 
        xDef = [ones(n,1) x]./sqrt(lambda(:,j,p));
        
        % Simulate observations from conditional posterior of beta and
        % sigma2 given lambda and the data.
        [estBeta(:,j + 1,p),estSigma2(1,j + 1,p)] = simulate(PriorMdl,xDef,yDef);
        
        % Estimate residuals.
        ep = y - [ones(n,1) x]*estBeta(:,j + 1,p);
        
        % Specify shape and scale using conditional posterior of lambda
        % given beta, sigma2, and the data
        sp = (nu(p) + 1)/2;
        sc =  2./(nu(p) + ep.^2/estSigma2(1,j + 1,p));
        
        % Draw from conditional posterior of lambda given beta, sigma2, 
        % and the data
        lambda(:,j + 1,p) = 1./gamrnd(sp,sc);
    end
end

Оцените следующие средние значения для $β$ , $σ^{2}$ , и $λ$ .

postEstBeta = squeeze(mean(estBeta(:,(burnin+1):end,:),2));
postLambda = squeeze(mean(lambda(:,(burnin+1):end,:),2));

Для каждого $ν$ , отобразите на графике данные и линии регрессии.

figure;
plot(x,y,'o');
h = gca;
xlim = h.XLim';
ylim = h.YLim;
hl = legend('Data');
hold on;
for p = 1:numNu;
    plotY = [ones(2,1) xlim]*postEstBeta(:,p);
    plot(xlim,plotY,'LineWidth',2);
    hl.String{p+1} = sprintf('nu = %f',nu(p));
end
xlabel('x');
ylabel('y');
title('Robust Bayesian Linear Regression')

Низкие значения $ν$ будьте склонны приписывать высокую изменчивость влиятельным выбросам. Поэтому линия регрессии напоминает строку regARIMA. Как $ν$ увеличения, строка ведет себя больше как те из наивного подхода.

Документация