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

Загрузите и предварительно обработайте данные

Загрузите США макроэкономический набор данных Data_USEconModel.mat.

load Data_USEconModel

Набор данных включает расписание MATLAB® DataTable, который содержит 14 переменных, измеренных от 1 квартала 1947 до 1 квартала 2009; UNRATE является уровнем безработицы США. Для получения дополнительной информации введите Description в командной строке.

Постройте весь ряд в той же фигуре, но в отдельных подграфиках.

figure;
for j = 1:size(DataTable,2)
    subplot(4,4,j)
    plot(DataTable.Time,DataTable{:,j});
    title(DataTable.Properties.VariableNames(j));
end

Все ряды кроме FEDFUNDS, GS10, TB3MS и UNRATE, кажется, имеют экспоненциальный тренд.

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

hasexpotrend = ~ismember(DataTable.Properties.VariableNames,...
    ["FEDFUNDS" "GD10" "TB3MS" "UNRATE"]);
DataTableLog = varfun(@log,DataTable,'InputVariables',...
    DataTable.Properties.VariableNames(hasexpotrend));
DataTableLog = [DataTableLog DataTable(:,DataTable.Properties.VariableNames(~hasexpotrend))];

DataTableLog является расписанием как DataTable, но те переменные с экспоненциальным трендом находятся на логарифмической шкале.

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

figure;
boxplot(DataTableLog.Variables,'Labels',DataTableLog.Properties.VariableNames);
h = gcf;
h.Position(3) = h.Position(3)*2.5;
title('Variable Box Plots');

Переменные имеют довольно подобные шкалы.

Чтобы определить, как хорошо модели временных рядов делают вывод, выполните эту процедуру:

Создайте оценку и предскажите демонстрационные переменные для данных о предикторе и ответа. Задайте горизонт прогноза 4 лет (16 четвертей).

fh = 16;
y = DataTableLog.UNRATE(1:(end - fh));
yF = DataTableLog.UNRATE((end - fh + 1):end);
isresponse = DataTable.Properties.VariableNames == "UNRATE";
X = DataTableLog{1:(end - fh),~isresponse};
XF = DataTableLog{(end - fh + 1):end,~isresponse};
p = size(X,2); % Number of predictors
predictornames = DataTableLog.Properties.VariableNames(~isresponse);

Соответствуйте простой модели линейной регрессии

Соответствуйте простой модели линейной регрессии к выборке оценки. Задайте имена переменных (имя переменной отклика должно быть последним элементом). Извлеките p-значения. Идентифицируйте незначительные коэффициенты при помощи 5%-го уровня значения.

SLRMdlFull = fitlm(X,y,'VarNames',DataTableLog.Properties.VariableNames)
slrPValue = SLRMdlFull.Coefficients.pValue;
isInSig = SLRMdlFull.CoefficientNames(slrPValue > 0.05)

SLRMdlFull = 


Linear regression model:
    UNRATE ~ [Linear formula with 14 terms in 13 predictors]

Estimated Coefficients:
                    Estimate       SE        tStat       pValue  
                    ________    ________    _______    __________

    (Intercept)       88.014      5.3229     16.535    2.6568e-37
    log_COE           7.1111      2.3426     3.0356     0.0027764
    log_CPIAUCSL     -3.4032      2.4611    -1.3828       0.16853
    log_GCE            -5.63      1.1581    -4.8615    2.6245e-06
    log_GDP           14.522      3.9406     3.6851    0.00030659
    log_GDPDEF        11.926      2.9298     4.0704    7.1598e-05
    log_GPDI         -2.5939     0.54603    -4.7504    4.2792e-06
    log_GS10         0.57467     0.29313     1.9605      0.051565
    log_HOANBS       -32.861      1.4507    -22.652    2.3116e-53
    log_M1SL         -3.3023     0.53185     -6.209     3.914e-09
    log_M2SL         -1.3845      1.0438    -1.3264       0.18647
    log_PCEC          -7.143      3.4302    -2.0824      0.038799
    FEDFUNDS        0.044542    0.047747    0.93287        0.3522
    TB3MS            -0.1577    0.064421     -2.448      0.015375


Number of observations: 185, Error degrees of freedom: 171
Root Mean Squared Error: 0.312
R-squared: 0.957,  Adjusted R-Squared: 0.954
F-statistic vs. constant model: 292, p-value = 2.26e-109

isInSig =

  1x4 cell array

    {'log_CPIAUCSL'}    {'log_GS10'}    {'log_M2SL'}    {'FEDFUNDS'}

SLRMdlFull является объектом модели LinearModel. P-значения предполагают, что на 5%-м уровне значения CPIAUCSL, GS10, M2SL и FEDFUNDS являются незначительными переменными.

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

idxreduced = ~ismember(predictornames,isInSig);
XReduced = X(:,idxreduced);
SLRMdlReduced = fitlm(XReduced,y,'VarNames',[predictornames(idxreduced) {'UNRATE'}]);

yFSLRF = predict(SLRMdlFull,XF);
yFSLRR = predict(SLRMdlReduced,XF(:,idxreduced));

fmseSLRF = sqrt(mean((yF - yFSLRF).^2))
fmseSLRR = sqrt(mean((yF - yFSLRR).^2))

fmseSLRF =

    0.6674


fmseSLRR =

    0.6105

FMSE упрощенной модели является меньше, чем FMSE полной модели, которая предполагает, что упрощенная модель делает вывод лучше.

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

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

[LassoBetaEstimates,FitInfo] = lasso(X,y,'Standardize',false,...
    'PredictorNames',predictornames);

По умолчанию lasso соответствует временам модели 100 регрессии лассо при помощи последовательности значений для. Поэтому LassoBetaEstimates 13 100 матрица оценок коэффициента регрессии. Строки соответствуют переменным прогноза в X, и столбцы соответствуют значениям. FitInfo является структурой, которая включает значения (FitInfo.Lambda), среднеквадратической ошибки для каждой подгонки (FitInfo.MSE) и предполагаемые прерывания (FitInfo.Intercept).

Вычислите FMSE каждой модели, возвращенной lasso.

yFLasso = FitInfo.Intercept + XF*LassoBetaEstimates;
fmseLasso = sqrt(mean((yF - yFLasso).^2,1));

Постройте значение коэффициентов регрессии относительно значения уменьшения.

figure;
hax = lassoPlot(LassoBetaEstimates,FitInfo);
L1Vals = hax.Children.XData;
yyaxis right
h = plot(L1Vals,fmseLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Frequentist Lasso')

Модель с 7 предикторами (df = 7), кажется, балансирует минимальный FMSE и сложность модели хорошо. Получите коэффициенты, соответствующие модели, содержащей 7 предикторов и приводящей к минимальному FMSE.

fmsebestlasso = min(fmseLasso(FitInfo.DF == 7));
idx = fmseLasso == fmsebestlasso;
bestLasso = [FitInfo.Intercept(idx); LassoBetaEstimates(:,idx)];
table(bestLasso,'RowNames',["Intercept" predictornames])

ans =

  14x1 table

                    bestLasso
                    _________

    Intercept         61.154 
    log_COE         0.042061 
    log_CPIAUCSL           0 
    log_GCE                0 
    log_GDP                0 
    log_GDPDEF         8.524 
    log_GPDI               0 
    log_GS10          1.6957 
    log_HOANBS       -18.937 
    log_M1SL         -1.3365 
    log_M2SL         0.17948 
    log_PCEC               0 
    FEDFUNDS               0 
    TB3MS           -0.14459 

Частотный анализ лассо предполагает, что переменные CPIAUCSL, GCE, GDP, GPDI, PCEC и FEDFUNDS или незначительны или избыточны.

Соответствуйте байесовой модели регрессии лассо

В представлении Bayesian регрессии лассо предшествующим распределением коэффициентов регрессии является Лаплас (удвойте экспоненциал), со средним значением 0 и шкалой, где фиксированный параметр уменьшения и. Для получения дополнительной информации смотрите lassoblm.

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

Когда вы реализуете Байесовую регрессию лассо в MATLAB®, знать о нескольких различиях между Statistics and Machine Learning Toolbox™ функционируют lasso и объект Econometrics Toolbox™ lassoblm и его присоединенные функции.

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

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

PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',predictornames);
table(PriorMdl.Lambda,'RowNames',PriorMdl.VarNames)

ans =

  14x1 table

                    Var1
                    ____

    Intercept       0.01
    log_COE            1
    log_CPIAUCSL       1
    log_GCE            1
    log_GDP            1
    log_GDPDEF         1
    log_GPDI           1
    log_GS10           1
    log_HOANBS         1
    log_M1SL           1
    log_M2SL           1
    log_PCEC           1
    FEDFUNDS           1
    TB3MS              1

PriorMdl является объектом модели lassoblm. lassoblm приписывает уменьшение 1 для каждого коэффициента кроме прерывания, которое имеет уменьшение 0.01. Эти значения по умолчанию не могут быть полезными в большинстве приложений; лучшая практика состоит в том, чтобы экспериментировать со многими значениями.

Считайте путь к регуляризации возвращенным lasso. Раздуйте значения уменьшения фактором, куда MSE лассо запущен, = 1 через количество значений уменьшения.

ismissing = any(isnan(X),2);
n = sum(~ismissing); % Effective sample size
lambda = FitInfo.Lambda*n./sqrt(FitInfo.MSE);

Считайте путь к регуляризации возвращенным lasso. Цикл через значения уменьшения в каждой итерации:

numlambda = numel(lambda);

% Preallocate
BayesLassoCoefficients = zeros(p+1,numlambda);
BayesLassoCI95 = zeros(p+1,2,numlambda);
fmseBayesLasso = zeros(numlambda,1);
BLCPlot = zeros(p+1,numlambda);

% Estimate and forecast
rng(10); % For reproducibility
for j = 1:numlambda
    PriorMdl.Lambda = lambda(j);
    [EstMdl,Summary] = estimate(PriorMdl,X,y,'Display',false);
    BayesLassoCoefficients(:,j) = Summary.Mean(1:(end - 1));
    BLCPlot(:,j) = Summary.Mean(1:(end - 1));
    BayesLassoCI95(:,:,j) = Summary.CI95(1:(end - 1),:);
    idx = BayesLassoCI95(:,2,j) > 0 & BayesLassoCI95(:,1,j) <= 0;
    BLCPlot(idx,j) = 0;
    yFBayesLasso = forecast(EstMdl,XF);
    fmseBayesLasso(j) = sqrt(mean((yF - yFBayesLasso).^2));
end

В каждой итерации estimate запускает сэмплер Гиббса к выборке от полных условных выражений (см. Аналитически Послушное Последующее поколение), и возвращает модель empiricalblm EstMdl, который содержит ничьи и сводную таблицу оценки Summary. Можно задать параметры для сэмплера Гиббса, такие как количество ничьих или утончающейся схемы.

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

Постройте следующие средние значения коэффициентов относительно нормированных штрафов L1 (сумма значений коэффициентов кроме прерывания). На том же графике, но в правильной оси y, строят FMSEs.

L1Vals = sum(abs(BLCPlot(2:end,:)),1)/max(sum(abs(BLCPlot(2:end,:)),1));

figure;
plot(L1Vals,BLCPlot(2:end,:))
xlabel('L1');
ylabel('Coefficient Estimates');
yyaxis right
h = plot(L1Vals,fmseBayesLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Bayesian Lasso')

Модель имеет тенденцию делать вывод лучше как нормированные увеличения штрафа L1 прошлые 0.1. Чтобы сбалансировать минимальный FMSE и сложность модели, выберите самую простую модель с FMSE, самым близким к 0,68.

idx = find(fmseBayesLasso > 0.68,1);
fmsebestbayeslasso = fmseBayesLasso(idx);
bestBayesLasso = BayesLassoCoefficients(:,idx);
bestCI95 = BayesLassoCI95(:,:,idx);
table(bestBayesLasso,bestCI95,'RowNames',["Intercept" predictornames])

ans =

  14x2 table

                    bestBayesLasso           bestCI95       
                    ______________    ______________________

    Intercept            90.587          79.819       101.62
    log_COE              6.5591          1.8093       11.239
    log_CPIAUCSL        -2.2971         -6.9019       1.7057
    log_GCE             -5.1707         -7.4902      -2.8583
    log_GDP              9.8839          2.3848       17.672
    log_GDPDEF           10.281          5.1677       15.936
    log_GPDI            -2.0973         -3.2168     -0.96728
    log_GS10            0.82485         0.22748       1.4237
    log_HOANBS          -32.538         -35.589      -29.526
    log_M1SL            -3.0656         -4.1499      -2.0066
    log_M2SL            -1.1007         -3.1243      0.75844
    log_PCEC            -3.3342         -9.6496       1.8482
    FEDFUNDS           0.043672       -0.056192      0.14254
    TB3MS              -0.15682        -0.29385    -0.022145

Один способ определить, незначительна ли переменная или избыточная, путем проверки, является ли 0 в его соответствующем коэффициенте 95%-м вероятным интервалом. Используя этот критерий, CPIAUCSL, M2SL, PCEC и FEDFUNDS незначительны.

Отобразите все оценки в таблице для сравнения.

SLRR = zeros(p + 1,1);
SLRR([true idxreduced]) = SLRMdlReduced.Coefficients.Estimate;
table([SLRMdlFull.Coefficients.Estimate; fmseSLRR],...
    [SLRR; fmseSLRR],...
    [bestLasso; fmsebestlasso],...
    [bestBayesLasso; fmsebestbayeslasso],'VariableNames',...
    {'SLRFull' 'SLRReduced' 'Lasso' 'BayesLasso'},...
    'RowNames',[PriorMdl.VarNames; 'MSE'])

ans =

  15x4 table

                    SLRFull     SLRReduced     Lasso      BayesLasso
                    ________    __________    ________    __________

    Intercept         88.014       88.327       61.154       90.587 
    log_COE           7.1111       10.854     0.042061       6.5591 
    log_CPIAUCSL     -3.4032            0            0      -2.2971 
    log_GCE            -5.63      -6.1958            0      -5.1707 
    log_GDP           14.522       16.701            0       9.8839 
    log_GDPDEF        11.926       9.1106        8.524       10.281 
    log_GPDI         -2.5939      -2.6963            0      -2.0973 
    log_GS10         0.57467            0       1.6957      0.82485 
    log_HOANBS       -32.861      -33.782      -18.937      -32.538 
    log_M1SL         -3.3023      -2.8099      -1.3365      -3.0656 
    log_M2SL         -1.3845            0      0.17948      -1.1007 
    log_PCEC          -7.143      -14.207            0      -3.3342 
    FEDFUNDS        0.044542            0            0     0.043672 
    TB3MS            -0.1577    -0.078449     -0.14459     -0.15682 
    MSE              0.61048      0.61048      0.79425      0.69639 

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

bestLambda = lambda(idx);
PriorMdl = bayeslm(p,'ModelType','lasso','Lambda',bestLambda,...
    'VarNames',predictornames);
PosteriorMdl = estimate(PriorMdl,[X; XF],[y; yF]);

Method: lasso MCMC sampling with 10000 draws
Number of observations: 201
Number of predictors:   14
 
              |   Mean      Std          CI95         Positive  Distribution 
-----------------------------------------------------------------------------
 Intercept    |  85.9643  5.2710   [75.544, 96.407]     1.000     Empirical  
 log_COE      |  12.3574  2.1817   [ 8.001, 16.651]     1.000     Empirical  
 log_CPIAUCSL |   0.1498  1.9150   [-3.733,  4.005]     0.525     Empirical  
 log_GCE      |  -7.1850  1.0556   [-9.208, -5.101]     0.000     Empirical  
 log_GDP      |  11.8648  3.9955   [ 4.111, 19.640]     0.999     Empirical  
 log_GDPDEF   |   8.8777  2.4221   [ 4.033, 13.745]     1.000     Empirical  
 log_GPDI     |  -2.5884  0.5294   [-3.628, -1.553]     0.000     Empirical  
 log_GS10     |   0.6910  0.2577   [ 0.194,  1.197]     0.997     Empirical  
 log_HOANBS   | -32.3356  1.4941  [-35.276, -29.433]    0.000     Empirical  
 log_M1SL     |  -2.2279  0.5043   [-3.202, -1.238]     0.000     Empirical  
 log_M2SL     |   0.3617  0.9123   [-1.438,  2.179]     0.652     Empirical  
 log_PCEC     | -11.3018  3.0353  [-17.318, -5.252]     0.000     Empirical  
 FEDFUNDS     |  -0.0132  0.0490   [-0.109,  0.082]     0.392     Empirical  
 TB3MS        |  -0.1128  0.0660   [-0.244,  0.016]     0.042     Empirical  
 Sigma2       |   0.1186  0.0122   [ 0.097,  0.145]     1.000     Empirical