Симулируйте данные

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

где:

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.2246e-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 распределил отклонение воздействия профилируемым параметром степеней свободы $$\nu$$Пусть:

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

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

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

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

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

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

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

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

Оцените следующие средние значения для $$\beta$$, $${\sigma }^2$$, и $$\lambda$$.

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

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

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')

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