Моделирование данных

Сгенерируйте случайную выборку из 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 итераций и примените период горения 5000. Предварительно выделите апостериорные рисунки и инициализируйте $$\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$$ увеличивается, линия ведет себя больше как линии наивного подхода.