simulate

Симулируйте коэффициенты регрессии и отклонение воздействия Байесовой модели линейной регрессии

Описание

пример

[BetaSim,sigma2Sim] = simulate(Mdl) возвращает случайный вектор из коэффициентов регрессии (BetaSim) и случайное отклонение воздействия (sigma2Sim) чертивший из Байесовой модели Mdl линейной регрессии из β и σ2.

  • Если Mdl объединенная предшествующая модель (возвращенный bayeslmто simulate чертит от предшествующих распределений.

  • Если Mdl объединенная следующая модель (возвращенный estimateто simulate чертит от апостериорных распределений.

пример

[BetaSim,sigma2Sim] = simulate(Mdl,X,y) чертит от крайних апостериорных распределений, произведенных или обновленных путем слияния данных о предикторе X и соответствующие данные об ответе y.

  • Если Mdl объединенная предшествующая модель, затем simulate производит крайние апостериорные распределения путем обновления предшествующей модели с информацией о параметрах, которые это получает из данных.

  • Если Mdl крайняя следующая модель, затем simulate обновляет последующее поколение с информацией о параметрах, которые она получает из дополнительных данных. Вероятность полных данных состоит из дополнительных данных X и y, и данные, которые создали Mdl.

NaNs в данных указывают на отсутствующие значения, который simulate удаляет при помощи мудрого списком удаления.

пример

[BetaSim,sigma2Sim] = simulate(___,Name,Value) использование любая из комбинаций входных аргументов в предыдущих синтаксисах и дополнительных опциях задано одним или несколькими аргументами пары "имя-значение". Например, можно задать значение для β или σ2 симулировать от апостериорного распределения conditional одного параметра, учитывая заданное значение другого параметра.

пример

[BetaSim,sigma2Sim,RegimeSim] = simulate(___) также возвращается, чертит от скрытого распределения режима если Mdl Байесова модель линейной регрессии для стохастического поискового выбора переменной (SSVS), то есть, если Mdl mixconjugateblm или mixsemiconjugateblm объект модели.

Примеры

свернуть все

Рассмотрите модель многофакторной линейной регрессии, которая предсказывает США действительный валовой национальный продукт (GNPR) при помощи линейной комбинации индекса промышленного производства (IPI), общая занятость (E), и действительная заработная плата (WR).

GNPRt=β0+β1IPIt+β2Et+β3WRt+εt.

\forall t, εt серия независимых Гауссовых воздействий со средним значением 0 и отклонение σ2.

Примите эти предшествующие распределения:

  • β|σ2N4(M,σ2V). M 4 1 вектор из средних значений, и V масштабированная положительная определенная ковариационная матрица 4 на 4.

  • σ2IG(A,B). A и B форма и шкала, соответственно, обратного гамма распределения.

Эти предположения и вероятность данных подразумевают нормальную обратную гамму сопряженная модель.

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.

load Data_NelsonPlosser
varNames = {'IPI' 'E' 'WR'};
X = DataTable{:,varNames};
y = DataTable{:,'GNPR'};

Создайте сопряженную предшествующую модель нормальной обратной гаммы для параметров линейной регрессии. Задайте количество предикторов p и имена переменных.

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',varNames);

PriorMdl conjugateblm Байесов объект модели линейной регрессии, представляющий предшествующее распределение коэффициентов регрессии и отклонения воздействия.

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

rng(1); % For reproducibility
[betaSimPrior,sigma2SimPrior] = simulate(PriorMdl)
betaSimPrior = 4×1

  -33.5917
  -49.1445
  -37.4492
  -25.3632

sigma2SimPrior = 0.1962

betaSimPrior случайным образом чертивший 4 1 вектор из коэффициентов регрессии, соответствующих именам в PriorMdl.VarNames. sigma2SimPrior выход является случайным образом чертившим скалярным отклонением воздействия.

Оцените апостериорное распределение.

PosteriorMdl = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood: -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

PosteriorMdl conjugateblm Байесов объект модели линейной регрессии, представляющий апостериорное распределение коэффициентов регрессии и отклонения воздействия.

Симулируйте набор коэффициентов регрессии и значение отклонения воздействия от апостериорного распределения.

[betaSimPost,sigma2SimPost] = simulate(PosteriorMdl)
betaSimPost = 4×1

  -25.9351
    4.4379
    0.0012
    2.4072

sigma2SimPost = 41.9575

betaSimPost и sigma2SimPost имейте те же размерности как betaSimPrior и sigma2SimPrior, соответственно, но чертятся от следующего.

Полагайте, что модель регрессии в Симулирует Значение параметров от Предшествующих и Апостериорных распределений.

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

load Data_NelsonPlosser
varNames = {'IPI' 'E' 'WR'};
X = DataTable{:,varNames};
y = DataTable{:,'GNPR'};
p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',varNames);
[PosteriorMdl,Summary] = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood: -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

Summary таблица, содержащая статистику что estimate отображения в командной строке.

Несмотря на то, что крайние и условные апостериорные распределения β и σ2 аналитически послушны, этот пример особое внимание о том, как реализовать сэмплер Гиббса, чтобы воспроизвести известные результаты.

Оцените модель снова, но используйте сэмплер Гиббса. Альтернатива между выборкой от условных апостериорных распределений параметров. Демонстрационные 10,000 раз и создайте переменные для предварительного выделения. Запустите сэмплер путем рисования от условного выражения, следующего из β данный σ2=2.

m = 1e4;
BetaDraws = zeros(p + 1,m);
sigma2Draws = zeros(1,m + 1);
sigma2Draws(1) = 2;
rng(1); % For reproducibility
for j = 1:m
    BetaDraws(:,j) = simulate(PriorMdl,X,y,'Sigma2',sigma2Draws(j));
    [~,sigma2Draws(j + 1)] = simulate(PriorMdl,X,y,'Beta',BetaDraws(:,j));
end
sigma2Draws = sigma2Draws(2:end); % Remove initial value from MCMC sample

Графики трассировки графика параметров.

figure;
for j = 1:(p + 1);
    subplot(2,2,j);
    plot(BetaDraws(j,:))
    ylabel('MCMC Draw')
    xlabel('Simulation Index')
    title(sprintf('Trace Plot — %s',PriorMdl.VarNames{j}));
end

Figure contains 4 axes objects. Axes object 1 with title Trace Plot — Intercept contains an object of type line. Axes object 2 with title Trace Plot — IPI contains an object of type line. Axes object 3 with title Trace Plot — E contains an object of type line. Axes object 4 with title Trace Plot — WR contains an object of type line.

figure;
plot(sigma2Draws)
ylabel('MCMC Draw')
xlabel('Simulation Index')
title('Trace plot — Sigma2')

Figure contains an axes object. The axes object with title Trace plot — Sigma2 contains an object of type line.

Выборки Цепи Маркова Монте-Карло (MCMC), кажется, сходятся и смешиваются хорошо.

Применяйтесь электротермотренировка 1 000 чертит, и затем вычислите средние значения и стандартные отклонения выборок MCMC. Сравните их с оценками от estimate.

bp = 1000;
postBetaMean = mean(BetaDraws(:,(bp + 1):end),2);
postSigma2Mean = mean(sigma2Draws(:,(bp + 1):end));
postBetaStd = std(BetaDraws(:,(bp + 1):end),[],2);
postSigma2Std = std(sigma2Draws((bp + 1):end));
[Summary(:,1:2),table([postBetaMean; postSigma2Mean],...
    [postBetaStd; postSigma2Std],'VariableNames',{'GibbsMean','GibbsStd'})]
ans=5×4 table
                   Mean          Std        GibbsMean     GibbsStd 
                 _________    __________    _________    __________

    Intercept      -24.249        8.7821      -24.293         8.748
    IPI             4.3913        0.1414       4.3917       0.13941
    E            0.0011202    0.00032931    0.0011229    0.00032875
    WR              2.4683       0.34895       2.4654       0.34364
    Sigma2          44.135         7.802       44.011        7.7816

Оценки очень близки. Изменения MCMC составляют различия.

Полагайте, что модель регрессии в Симулирует Значение параметров от Предшествующих и Апостериорных распределений.

Примите эти предшествующие распределения для k = 0,...,3:

  • βk|σ2,γk=γkσVk1Z1+(1-γk)σVk2Z2, где Z1 и Z2 независимые, стандартные нормальные случайные переменные. Поэтому коэффициенты имеют Гауссово распределение смеси. Примите, что все коэффициенты условно независимы, априорно, но они зависят от отклонения воздействия.

  • σ2IG(A,B). A и B форма и шкала, соответственно, обратного гамма распределения.

  • γk{0,1}и это представляет переменную режима включения случайной переменной дискретным равномерным распределением.

Создайте предшествующую модель для выполнения SSVS. AssumeThat β и σ2зависят (сопряженная модель смеси). Задайте количество предикторов p и имена коэффициентов регрессии.

p = 3;
PriorMdl = mixconjugateblm(p,'VarNames',["IPI" "E" "WR"]);

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

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

cardRegime = 2^(PriorMdl.Intercept + PriorMdl.NumPredictors)
cardRegime = 16

Симулируйте 10 000 режимов от апостериорного распределения.

rng(1);
[~,~,RegimeSim] = simulate(PriorMdl,X,y,'NumDraws',10000);

RegimeSim 4 1000 логическая матрица. Строки соответствуют переменным в Mdl.VarNames, и столбцы соответствуют, чертит от апостериорного распределения.

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

cRegime = num2cell(RegimeSim,1);
cRegime = categorical(cellfun(@(c)join(PriorMdl.VarNames(c),"."),cRegime));
cRegime(ismissing(cRegime)) = "NoCoefficients";

histogram(cRegime);
title('Variables Included in Models')
ylabel('Frequency');

Figure contains an axes object. The axes object with title Variables Included in Models contains an object of type categoricalhistogram.

Вычислите крайнюю апостериорную вероятность переменного включения.

table(mean(RegimeSim,2),'RowNames',PriorMdl.VarNames,...
    'VariableNames',"Regime")
ans=4×1 table
                 Regime
                 ______

    Intercept    0.8829
    IPI          0.4547
    E             0.098
    WR           0.1692

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

  • λjIG(ν/2,2/ν).

  • εj|λjN(0,λjσ2)

  • f(β,σ2)1σ2

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

  • εjt(0,σ2,ν)

  • λj|εjIG(ν+12,2ν+εj2/σ2)

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

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

Сгенерировать n=100 ответы от yt=1+2xt+et, где x[0,2] и etN(0,0.52).

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;

Подбор линейной модели к данным. Отобразите на графике данные и подходящую линию регрессии.

Mdl = fitlm(x,y)
Mdl = 
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
figure;
plot(Mdl);
hl = legend;
hold on;

Figure contains an axes object. The axes object with title y vs. x1 contains 4 objects of type line. These objects represent Data, Fit, Confidence bounds.

Симулированные выбросы, кажется, влияют на подходящую линию регрессии.

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

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

  2. Вычислите остаточные значения.

  3. Чертите значения от условного выражения, следующего из λ.

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

m = 20000;
nu = 1;
burnin = 5000;
lambda = ones(n,m + 1);
estBeta = zeros(2,m + 1);
estSigma2 = zeros(1,m + 1);
for j = 1:m
    yDef = y./sqrt(lambda(:,j));
    xDef = [ones(n,1) x]./sqrt(lambda(:,j));
    PriorMdl = bayeslm(2,'Model','diffuse','Intercept',false);
    [estBeta(:,j + 1),estSigma2(1,j + 1)] = simulate(PriorMdl,xDef,yDef);
    ep = y - [ones(n,1) x]*estBeta(:,j + 1);
    sp = (nu + 1)/2;
    sc =  2./(nu + ep.^2/estSigma2(1,j + 1));
    lambda(:,j + 1) = 1./gamrnd(sp,sc);
end

Хорошая практика должна диагностировать сэмплер MCMC путем исследования графиков трассировки. Для краткости этот пример пропускает эту задачу.

Вычислите среднее значение ничьих от следующих из коэффициентов регрессии. Удалите электротермотренировку, чертит.

postEstBeta = mean(estBeta(:,(burnin + 1):end),2)
postEstBeta = 2×1

    1.3971
    1.7051

Оценка точки пересечения ниже, и наклон выше, чем оценки, возвращенные fitlm.

Постройте устойчивую линию регрессии с линией регрессии, адаптированной наименьшими квадратами.

h = gca;
xlim = h.XLim';
plotY = [ones(2,1) xlim]*postEstBeta;
plot(xlim,plotY,'LineWidth',2);
hl.String{4} = 'Robust Bayes';

Figure contains an axes object. The axes object with title y vs. x1 contains 5 objects of type line. These objects represent Data, Fit, Confidence bounds, Robust Bayes.

Подгонка линии регрессии использование устойчивой Байесовой регрессии, кажется, лучшая подгонка.

Максимум по опыту вероятность (MAP), оценка является следующим режимом, то есть, значение параметров, которое дает к максимуму следующей PDF. Если следующее аналитически тяжело, то можно использовать выборку Монте-Карло, чтобы оценить MAP.

Полагайте, что модель линейной регрессии в Симулирует Значение параметров от Предшествующих и Апостериорных распределений.

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.

load Data_NelsonPlosser
varNames = {'IPI' 'E' 'WR'};
X = DataTable{:,varNames};
y = DataTable{:,'GNPR'};

Создайте сопряженную предшествующую модель нормальной обратной гаммы для параметров линейной регрессии. Задайте количество предикторов p и имена переменных.

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',varNames)
PriorMdl = 
  conjugateblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
               Mu: [4x1 double]
                V: [4x4 double]
                A: 3
                B: 1

 
           |  Mean     Std            CI95         Positive       Distribution     
-----------------------------------------------------------------------------------
 Intercept |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

Оцените крайние апостериорные распределения β и σ2.

rng(1); % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood: -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

Отображение включает крайнюю статистику апостериорного распределения.

Извлеките следующее среднее значение β из следующей модели и извлечения следующая ковариация β из сводных данных оценки, возвращенных summarize.

estBetaMean = PosteriorMdl.Mu;
Summary = summarize(PosteriorMdl);
EstBetaCov = Summary.Covariances{1:(end - 1),1:(end - 1)};

estBetaMean 4 1 вектор, представляющий среднее значение крайнего следующего из β. EstBetaCov матрица 4 на 4, представляющая ковариационную матрицу следующего из β.

Чертите 10 000 значений параметров от апостериорного распределения.

rng(1); % For reproducibility
[BetaSim,sigma2Sim] = simulate(PosteriorMdl,'NumDraws',1e5);

BetaSim 4 10,000 матрица случайным образом чертивших коэффициентов регрессии. sigma2Sim 1 10,000 вектор из случайным образом чертивших отклонений воздействия.

Транспонируйте и стандартизируйте матрицу коэффициентов регрессии. Вычислите корреляционную матрицу коэффициентов регрессии.

estBetaStd = sqrt(diag(EstBetaCov)');
BetaSim = BetaSim';
BetaSimStd = (BetaSim - estBetaMean')./estBetaStd;
BetaCorr = corrcov(EstBetaCov);
BetaCorr = (BetaCorr + BetaCorr')/2; % Enforce symmetry

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

betaPDF = mvtpdf(BetaSimStd,BetaCorr,68);
a = 34;
b = 0.00069;
igPDF = @(x,ap,bp)1./(gamma(ap).*bp.^ap).*x.^(-ap-1).*exp(-1./(x.*bp));...
    % Inverse gamma pdf
sigma2PDF = igPDF(sigma2Sim,a,b);

Найдите симулированные значения, которые максимизируют соответствующий pdfs, то есть, следующие режимы.

[~,idxMAPBeta] = max(betaPDF);
[~,idxMAPSigma2] = max(sigma2PDF);
betaMAP = BetaSim(idxMAPBeta,:);
sigma2MAP = sigma2Sim(idxMAPSigma2);

betaMAP и sigma2MAP оценки MAP.

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

table(betaMAP',PosteriorMdl.Mu,'VariableNames',{'MAP','Mean'},...
    'RowNames',PriorMdl.VarNames)
ans=4×2 table
                    MAP         Mean   
                 _________    _________

    Intercept      -24.559      -24.249
    IPI             4.3964       4.3913
    E            0.0011389    0.0011202
    WR              2.4473       2.4683

Оценки справедливо друг близко к другу.

Оцените аналитический режим следующего из σ2. Сравните его с предполагаемым MAP σ2.

igMode = 1/(b*(a+1))
igMode = 41.4079
sigma2MAP
sigma2MAP = 41.4075

Эти оценки также довольно близки.

Входные параметры

свернуть все

Стандартная Байесова модель линейной регрессии или модель для выбора переменного предиктора в виде объекта модели в этой таблице.

Объект моделиОписание
conjugateblmЗависимый, нормальная обратная гамма спрягает модель, возвращенную bayeslm или estimate
semiconjugateblmНезависимый, нормальная обратная гамма полуспрягает модель, возвращенную bayeslm
diffuseblmРассейте предшествующую модель, возвращенную bayeslm
empiricalblmПредшествующая модель, охарактеризованная выборками от предшествующих распределений, возвращенных bayeslm или estimate
customblmПредшествующая функция распределения, которую вы объявляете возвращенный bayeslm
mixconjugateblmЗависимая, Гауссова гамма инверсии смеси спрягает модель для выбора переменного предиктора SSVS, возвращенного bayeslm
mixsemiconjugateblmНезависимая, Гауссова гамма инверсии смеси полуспрягает модель для выбора переменного предиктора SSVS, возвращенного bayeslm
lassoblmБайесова модель регрессии лассо, возвращенная bayeslm

Примечание

  • Как правило, объекты модели, возвращенные estimate представляйте крайние апостериорные распределения. Когда вы оцениваете следующее при помощи estimate, если вы задаете оценку следующего условного выражения, то estimate возвращает предшествующую модель.

  • Если Mdl diffuseblm модель, затем необходимо также предоставить X и y потому что simulate не может чертить от неподходящего предшествующего распределения.

  • Если вы предоставляете lassoblm, mixconjugateblm, или mixsemiconjugateblm объект модели, снабдите данными X и y, и чертите одно значение от следующего, затем лучшая практика должна инициализировать сэмплер Гиббса путем определения BetaStart и Sigma2Start аргументы в виде пар имя-значение.

Данные о предикторе для модели многофакторной линейной регрессии в виде numObservations- PriorMdl.NumPredictors числовая матрица. numObservations количество наблюдений и должно быть равно длине y.

Если Mdl апостериорное распределение, затем столбцы X должен соответствовать столбцам данных о предикторе, используемых, чтобы оценить следующее.

Типы данных: double

Данные об ответе для модели многофакторной линейной регрессии в виде числового вектора с numObservations элементы.

Типы данных: double

Аргументы name-value

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: 'Sigma2',2 задает симуляцию от условного апостериорного распределения коэффициентов регрессии, учитывая данные и заданное отклонение воздействия 2.
Опции для всех моделей

свернуть все

Количество ничьих к выборке от распределения MdlВ виде разделенной запятой пары, состоящей из 'NumDraws' и положительное целое число.

Совет

Если Mdl empiricalblm или customblm объект модели, затем хорошая практика должна задать электротермотренировку с BurnIn и утончающийся множитель с Thin. Для получения дополнительной информации на настроенном объеме выборки, см. Алгоритмы.

Пример: 'NumDraws',1e7

Типы данных: double

Опции для всех моделей кроме эмпирического

свернуть все

Значение коэффициентов регрессии для симуляции от условного распределения отклонения воздействия в виде разделенной запятой пары, состоящей из 'Beta' и (Mdl.Intercept + Mdl.NumPredictors)-by-1 числовой вектор. При использовании апостериорного распределения, simulate чертит от π (σ2|yX, β = Beta), где y yX X, и Beta значение 'Beta'. Если Mdl.Intercept true, затем Beta(1) соответствует точке пересечения модели. Все другие значения соответствуют переменным предикторам, которые составляют столбцы X.

Вы не можете задать Beta и Sigma2 одновременно.

По умолчанию, simulate не чертит от условного выражения, следующего из σ2.

Пример: 'Beta',1:3

Типы данных: double

Значение отклонения воздействия для симуляции от условного распределения коэффициентов регрессии в виде разделенной запятой пары, состоящей из 'Sigma2' и положительный числовой скаляр. При использовании апостериорного распределения, simulate чертит от π (β |yX, Sigma2), где y yX X, и Sigma2 значение 'Sigma2'.

Вы не можете задать Sigma2 и Beta одновременно.

По умолчанию, simulate не чертит от условного выражения, следующего из β.

Пример: 'Sigma2',1

Типы данных: double

Опции для всех моделей кроме сопряженного и эмпирического

свернуть все

Количество ничьих, чтобы удалить с начала выборки уменьшать переходные эффекты в виде разделенной запятой пары, состоящей из 'BurnIn' и неотрицательный скаляр. Для получения дополнительной информации, на как simulate уменьшает полную выборку, см. Алгоритмы.

Совет

Помочь вам задать соответствующий размер электротермотренировки:

  1. Определите степень переходного поведения в выборке путем определения 'BurnIn',0.

  2. Симулируйте несколько тысяч наблюдений при помощи simulate.

  3. Постройте графики трассировки.

Пример: 'BurnIn',0

Типы данных: double

Настроенный множитель объема выборки в виде разделенной запятой пары, состоящей из 'Thin' и положительное целое число.

Фактическим объемом выборки является BurnIn + NumDraws*Thin. После отбрасывания выжигания дефектов, simulate отбрасывает каждый Thin– 1 чертит, и затем сохраняет следующую ничью. Для получения дополнительной информации о как simulate уменьшает полную выборку, см. Алгоритмы.

Совет

Уменьшать потенциальную большую последовательную корреляцию в выборке или уменьшать потребление памяти ничьих, сохраненных в PosteriorMdl, задайте большое значение для Thin.

Пример: 'Thin',5

Типы данных: double

Начальные значения коэффициентов регрессии для сэмплера в виде разделенной запятой пары, состоящей из 'BetaStart' и числовой вектор-столбец с (Mdl.Intercept + Mdl.NumPredictors) элементы. По умолчанию, BetaStart оценка обычных наименьших квадратов (OLS).

Совет

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

Пример: 'BetaStart',[1; 2; 3]

Типы данных: double

Начальные значения отклонения воздействия для сэмплера в виде разделенной запятой пары, состоящей из 'Sigma2Start' и положительный числовой скаляр. По умолчанию, Sigma2Start остаточная среднеквадратическая ошибка OLS.

Совет

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

Пример: 'Sigma2Start',4

Типы данных: double

Опции для моделей SSVS

свернуть все

Начальные значения скрытых режимов для сэмплера в виде разделенной запятой пары, состоящей из 'RegimeStart' и логический вектор-столбец с (Mdl.Intercept + Mdl.NumPredictors) элементы. RegimeStart (k) = true указывает на включение переменной Mdl. VarNames (k), и RegimeStart (k) = false указывает на исключение той переменной.

Совет

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

Пример: 'RegimeStart',logical(randi([0 1],Mdl.Intercept + Mdl.NumPredictors,1))

Типы данных: double

Опции для пользовательских моделей

свернуть все

Репараметризация σ2 как журнал (σ2) во время следующей оценки и симуляции в виде разделенной запятой пары, состоящей из 'Reparameterize' и значение в этой таблице.

ЗначениеОписание
falsesimulate не повторно параметрирует σ2.
truesimulate повторно параметрирует σ2 как журнал (σ2). simulate преобразует результаты назад в исходную шкалу и не изменяет функциональную форму PriorMdl.LogPDF.

Совет

Если вы испытываете числовую нестабильность во время следующей оценки или симуляции σ2, затем задайте 'Reparameterize',true.

Пример: 'Reparameterize',true

Типы данных: логический

Сэмплер MCMC в виде разделенной запятой пары, состоящей из 'Sampler' и значение в этой таблице.

ЗначениеОписание
'slice'Сэмплер среза
'metropolis'Случайный сэмплер Столицы обхода
'hmc'Сэмплер Гамильтонова Монте-Карло (HMC)

Совет

  • Увеличить качество MCMC чертит, настройте сэмплер.

    1. Перед вызовом simulate, задайте настраивающиеся параметры и их значения при помощи sampleroptions. Например, чтобы задать ширину сэмплера среза widthИспользование:

      options = sampleroptions('Sampler',"slice",'Width',width);

    2. Задайте объект, содержащий настраивающиеся технические требования параметра, возвращенные sampleroptions при помощи 'Options' аргумент пары "имя-значение". Например, чтобы использовать настраивающиеся технические требования параметра в options, задайте:

      'Options',options

  • Если вы задаете сэмплер HMC, то лучшая практика должна предоставить градиент для некоторых переменных, по крайней мере. simulate обращается численный расчет любых недостающих частных производных (NaN значения) в векторе градиента.

Пример: 'Sampler',"hmc"

Типы данных: string

Опции сэмплера в виде разделенной запятой пары, состоящей из 'Options' и массив структур, возвращенный sampleroptions. Используйте 'Options' задавать сэмплер MCMC и его настраивающие значения параметров.

Пример: 'Options',sampleroptions('Sampler',"hmc")

Типы данных: struct

Типичная ширина интервала выборки вокруг текущего значения в предельных распределениях для сэмплера среза в виде разделенной запятой пары, состоящей из 'Width' и положительный числовой скаляр или (PriorMdl.Intercept + PriorMdl.NumPredictors+ 1 )-by-1 числовой вектор из положительных значений. Первый элемент соответствует точке пересечения модели, если вы существуете в модели. Следующий PriorMdl.NumPredictors элементы соответствуют коэффициентам переменных предикторов, упорядоченных столбцами данных предиктора. Последний элемент соответствует отклонению модели.

  • Если Width скаляр, затем simulate применяет Width ко всему PriorMdl.NumPredictors + PriorMdl.Intercept+ 1 предельные распределения.

  • Если Width числовой вектор, затем simulate применяет первый элемент к точке пересечения (если вы существуете), следующий PriorMdl.NumPredictors элементы к коэффициентам регрессии, соответствующим переменным предикторам в X, и последний элемент к отклонению воздействия.

  • Если объем выборки (size(X,1)) меньше 100, затем Width 10 по умолчанию.

  • Если объем выборки - по крайней мере 100, то simulate наборы Width к вектору из соответствующих следующих стандартных отклонений по умолчанию, принимая рассеянную предшествующую модель (diffuseblm).

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

simulate отправляет Width к slicesample функция. Для получения дополнительной информации смотрите slicesample.

Совет

  • Для максимальной гибкости задайте ширину сэмплера среза width при помощи 'Options' аргумент пары "имя-значение". Например:

    'Options',sampleroptions('Sampler',"slice",'Width',width)

Пример: 'Width',[100*ones(3,1);10]

Выходные аргументы

свернуть все

Симулированные коэффициенты регрессии, возвращенные как (Mdl.Intercept + Mdl.NumPredictors)-by-NumDraws числовая матрица. Строки соответствуют переменным в Mdl.VarNames, и столбцы соответствуют отдельным, последовательным, независимым ничьим от распределения.

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

Симулированные режимы, указывающие на переменное включение или исключение из модели, возвращенной как (Mdl.Intercept + Mdl.NumPredictors)-by-NumDraws логическая матрица. Строки соответствуют переменным в Mdl.VarNames, и столбцы соответствуют отдельным, последовательным, независимым ничьим от распределения.

simulate возвращает Regime только если Mdl mixconjugateblm или mixsemiconjugateblm объект модели.

RegimeSim (kD) = true указывает на включение переменной Mdl. VarNames (k) в модели ничьей d, и RegimeSim (kD) = false указывает на исключение той переменной в модели ничьей d.

Ограничения

  • simulate не может чертить значения от improper distribution, то есть, распределение, плотность которого не объединяется к 1.

  • Если Mdl empiricalblm объект модели, затем вы не можете задать Beta или Sigma2. Вы не можете симулировать от условных апостериорных распределений при помощи эмпирического распределения.

Больше о

свернуть все

Байесова модель линейной регрессии

Bayesian linear regression model обрабатывает параметры β и σ2 в модели yt многофакторной линейной регрессии (MLR) = xt β + εt как случайные переменные.

В течение многих времен t = 1..., T:

  • yt является наблюдаемым ответом.

  • xt является 1 на (p + 1) вектор-строка из наблюдаемых величин предикторов p. Вмещать точку пересечения модели, x 1t = 1 для всего t.

  • β (p + 1)-by-1 вектор-столбец коэффициентов регрессии, соответствующих переменным, которые составляют столбцы xt.

  • εt является случайным воздействием со средним значением нуля и Cov (ε) = σ2T I ×T, в то время как ε является T-by-1 вектор, содержащий все воздействия. Эти предположения подразумевают, что вероятность данных

    (β,σ2|y,x)=t=1Tϕ(yt;xtβ,σ2).

    ϕ (yt; xtβ, σ2) Гауссова плотность вероятности со средним xtβ и отклонением σ2 оцененный в yt;.

Прежде, чем рассмотреть данные, вы налагаете предположение joint prior distribution на (β, σ2). В Байесовом анализе вы обновляете распределение параметров при помощи информации о параметрах, полученных из вероятности данных. Результатом является joint posterior distribution (β, σ2) или conditional posterior distributions параметров.

Алгоритмы

  • Каждый раз, когда simulate должен оценить апостериорное распределение (например, когда Mdl представляет предшествующее распределение, и вы предоставляете X и y) и следующее аналитически послушно, simulate симулирует непосредственно от следующего. В противном случае, simulate обращения к симуляции Монте-Карло, чтобы оценить следующее. Для получения дополнительной информации смотрите Следующую Оценку и Вывод.

  • Если Mdl объединенная следующая модель, затем simulate симулирует данные из него по-другому по сравнению с когда Mdl объединенная предшествующая модель, и вы предоставляете X и y. Поэтому, если вы устанавливаете тот же случайный seed и генерируете случайные значения оба пути, затем вы не можете получить те же значения. Однако соответствующие эмпирические распределения на основе достаточного числа ничьих эффективно эквивалентны.

  • Этот рисунок показывает как simulate уменьшает выборку при помощи значений NumDraws, Thin, и BurnIn.

    Прямоугольники представляют последовательные ничьи от распределения. simulate удаляет белые прямоугольники из выборки. Остающийся NumDraws черные прямоугольники составляют выборку.

  • Если Mdl semiconjugateblm объект модели, затем simulate выборки от апостериорного распределения путем применения сэмплера Гиббса.

    1. simulate использует значение по умолчанию Sigma2Start для σ2 и чертит значение β от π (β |σ2x, y .

    2. simulate чертит значение σ2 от π (σ2|β, X, y) при помощи ранее сгенерированного значения β.

    3. Функция повторяет шаги 1 и 2 до сходимости. Чтобы оценить сходимость, постройте график трассировки выборки.

    Если вы задаете BetaStartто simulate чертит значение σ2 от π (σ2|β, X, y), чтобы запустить сэмплер Гиббса. simulate не возвращает это сгенерированное значение σ2.

  • Если Mdl empiricalblm объект модели и вы не предоставляете X и yто simulate чертит от Mdl.BetaDraws и Mdl.Sigma2Draws. Если NumDraws меньше чем или равно numel(Mdl.Sigma2Draws)то simulate возвращает первый NumDraws элементы Mdl.BetaDraws и Mdl.Sigma2Draws как случайные ничьи для соответствующего параметра. В противном случае, simulate случайным образом передискретизирует NumDraws элементы от Mdl.BetaDraws и Mdl.Sigma2Draws.

  • Если Mdl customblm объект модели, затем simulate использует сэмплер MCMC, чтобы чертить от апостериорного распределения. В каждой итерации программное обеспечение конкатенирует текущие значения коэффициентов регрессии и отклонения воздействия в (Mdl.Intercept + Mdl.NumPredictors + 1) вектор-by-1 и передачи это к Mdl.LogPDF. Значение отклонения воздействия является последним элементом этого вектора.

  • Сэмплер HMC требует и логарифмической плотности и ее градиента. Градиентом должен быть (NumPredictors+Intercept+1)- 1 вектор. Если производные определенных параметров затрудняют, чтобы вычислить, то, в соответствующих местоположениях градиента, NaN предоставления значения вместо этого. simulate замены NaN значения с числовыми производными.

  • Если Mdl lassoblm, mixconjugateblm, или mixsemiconjugateblm объект модели и вы предоставляете X и yто simulate выборки от апостериорного распределения путем применения сэмплера Гиббса. Если вы не снабжаете данными, то simulate выборки от аналитических, безусловных предшествующих распределений.

  • simulate не возвращает начальные значения по умолчанию, которые это генерирует.

  • Если Mdl mixconjugateblm или mixsemiconjugateblmто simulate чертит от распределения режима сначала, учитывая текущее состояние цепи (значения RegimeStart, BetaStart, и Sigma2Start). Если вы чертите одну выборку и не задаете значения для RegimeStart, BetaStart, и Sigma2Startто simulate использует значения по умолчанию и выдает предупреждение.

Введенный в R2017a