Estimador tendencioso para regressão, obtendo melhores resultados do que imparcial no modelo Erro nas variáveis

Estou trabalhando em alguns dados sintáticos para o modelo Error In Variable para algumas pesquisas. Atualmente, tenho uma única variável independente e estou assumindo que conheço a variação do valor verdadeiro da variável dependente.

Portanto, com essas informações, posso obter um estimador imparcial para o coeficiente da variável dependente.

O modelo:

$\tilde{x} = x + e_1$
$y = 0.5x -10 + e_2$
Onde:
$e_1\text{~}N(0,\sigma^2)$ para alguns $\sigma$
$e_2\text{~}N(0,1)$

Onde os valores de $y,\tilde{x}$ são conhecidos apenas para cada amostra e também o desvio padrão do valor real de $x$ para a amostra: $\sigma_x$ .

Recebo o tendenciosa coeficiente) usando OLS, e depois fazer ajustes usando: $\hat{\beta}$

$\beta' = \hat{\beta} * \frac{\hat{\sigma}_\tilde{x}^2}{\sigma_x^2}$

Vejo que meu novo estimador imparcial para o coeficiente é muito melhor (mais próximo do valor real) com esse modelo, mas o MSE está ficando pior do que usar o estimador enviesado.

O que está acontecendo? Eu esperava que um estimador tendencioso produzisse melhores resultados que o tendencioso.

Código Matlab:

reg_mse_agg = [];
fixed_mse_agg = [];
varMult = 1;
numTests = 60;
for dataNumber=1:8
    reg_mses = [];
    fixed_mses = [];

    X = rand(1000,1);
    X(:,1) = X(:,1) * 10;
    X(:,1) = X(:,1) + 5;

    varX = var(X);
    y = 0.5 * X(:,1) -10;

    y = y + normrnd(0,1,size(y));    
    origX = X;
    X = X + normrnd(0,dataNumber * varMult ,size(X));
    train_size = floor(0.5 * length(y));
    for t=1:numTests,
        idx = randperm(length(y));
        train_idx = idx(1:train_size);
        test_idx = idx(train_size+1:end);
        Xtrain = X(train_idx,:);
        ytrain = y(train_idx);        
        Xtest = X(test_idx,:);
        ytest = y(test_idx);

        b = OLS_solver(Xtrain, ytrain);
        %first arg of evaluate returns MSE, working correctly.
        [ reg_mse, ~ ] = evaluate( b,Xtest,ytest);
        reg_mses = [reg_mses ; reg_mse];

        varInd = var(Xtrain);
        varNoise = varInd - varX;
        bFixed = [0 0]';
        bFixed(1) = b(1) * varInd / varX;
        bFixed(2) = mean(ytrain - bFixed(1)*Xtrain);
        [fixed_mse,~ ] = evaluate( bFixed,Xtest,ytest);
        fixed_mses = [fixed_mses ; fixed_mse];

        dataNumber * varMult
        b
        bFixed

    end
    reg_mse_agg = [reg_mse_agg , reg_mses];
    fixed_mse_agg = [fixed_mse_agg , fixed_mses];

end

mean(reg_mse_agg)
mean(fixed_mse_agg)

Resultados:

MSE do estimador enviesado:

ans =

  Columns 1 through 7

    1.2171    1.6513    1.9989    2.3914    2.5766    2.6712    2.5997

  Column 8

    2.8346

MSE do estimador imparcial:

ans =

  Columns 1 through 7

    1.2308    2.0001    2.9555    4.9727    7.6757   11.3106   14.4283

  Column 8

   11.5653

Além disso, imprimir os valores de be bFixed- vejo que bFixedé realmente mais próximo dos valores reais do 0.5,-10que o estimador enviesado (como esperado).

PS Os resultados da imparcialidade ser pior do que o estimador enviesado são estatisticamente significativos - o teste é omitido no código, pois é uma simplificação do código da "versão completa".

UPDTAE: Eu adicionado um teste que verifica e , e o estimador parcial realmente é significativamente pior (valor maior) do que a um imparcial de acordo com a essa métrica, mesmo que o MSE do estimador enviesado (no conjunto de testes) seja significativamente melhor. Onde é o coeficiente real da variável é o estimador enviesada para , e é o estimador de $\sum_{\text{for each test}}{(\hat{\beta}-\beta)^2}$ $\sum_{\text{for each test}}{(\beta'-\beta)^2}$
$\beta=0.5$ $\hat{\beta}$ $\beta$ $\beta'$ $\beta$ .

Acredito que isso mostre que a razão para os resultados NÃO é a maior variação do estimador imparcial, pois ainda está mais próximo do valor real.

Crédito: Usando notas da palestra de Steve Pischke como recurso

regression matlab unbiased-estimator errors-in-variables amit
fonte

Seria útil se você postasse também os resultados, e não apenas o código.

Alecos Papadopoulos 31/03

@AlecosPapadopoulos Adicionou, não adicionou a impressão de todos os valores de be bFixed, mas explicou o que eles mostram.

Am 31/03

Respostas:

Resumo: os parâmetros corrigidos são para previsão em função do verdadeiro preditor . Se for usado na previsão, os parâmetros originais terão melhor desempenho. $x$ $\tilde{x}$

Observe que existem dois modelos de previsão linear diferentes à espreita. Em primeiro lugar, dada , $y$ $x$ segundo, dado ,

{\hat{y}}_{x} = β x + α,

$\begin{equation} \hat{y}_x = \beta\,x + \alpha, \end{equation}$

y

$y$

\tilde{x}

$\tilde{x}$

{\hat{y}}_{\tilde{x}} = \tilde{β} \tilde{x} + \tilde{α} .

$\begin{equation} \hat{y}_{\tilde{x}} = \tilde{\beta}\,\tilde{x} + \tilde{\alpha}. \end{equation}$

$x$ $\tilde{x}$

$\hat{\tilde{\beta}},\hat{\tilde{\alpha}}$
$\hat{\beta},\hat{\alpha}$
$\hat{y}_1 = \hat{\beta}\,\tilde{x} + \hat{\alpha}$ $\hat{y}_2 = \hat{\tilde{\beta}}\,\tilde{x} + \hat{\tilde{\alpha}}$

$\tilde{x}$ $x$

$\alpha$ $\beta$ $\tilde{x}$ $x$ $x$

\hat{\hat{y_{x}}} = β \hat{x} (\tilde{x}) + α = β (μ_{x} + (\hat{x} - μ_{x}) \frac{σ_{x}^{2}}{σ_{\tilde{x}}^{2}}) + α = \frac{σ_{x}^{2}}{σ_{\tilde{x}}^{2}} β + α - β (1 - \frac{σ_{x}^{2}}{σ_{\tilde{x}}^{2}}) μ_{x} .

$\begin{equation} \hat{\hat{y_x}} = \beta\,\hat{x}(\tilde{x}) + \alpha = \beta\, (\mu_x + (\hat{x}-\mu_x)\,\frac{\sigma^2_{x}}{\sigma^2_{\tilde{x}}}) + \alpha = \frac{\sigma_x^2}{\sigma^2_{\tilde{x}}}\beta + \alpha - \beta(1-\frac{\sigma_x^2}{\sigma^2_{\tilde{x}}})\mu_x. \end{equation}$

\tilde{α}, \tilde{β}

$\tilde{\alpha},\tilde{\beta}$

α, β

$\alpha,\beta$

\tilde{α}, \tilde{β}

$\tilde{\alpha},\tilde{\beta}$

α, β

$\alpha,\beta$

Teste

Editei o código no OP para avaliar também as MPEs para previsões usando a versão não barulhenta da previsão (código no final da resposta). Os resultados são

Reg parameters, noisy predictor
1.3387    1.6696    2.1265    2.4806    2.5679    2.5062    2.5160    2.8684

Fixed parameters, noisy predictor
1.3981    2.0626    3.2971    5.0220    7.6490   10.2568   14.1139   20.7604

Reg parameters, true predictor
1.3354    1.6657    2.1329    2.4885    2.5688    2.5198    2.5085    2.8676

Fixed parameters, true predictor
1.1139    1.0078    1.0499    1.0212    1.0492    0.9925    1.0217    1.2528

$x$ $\tilde{x}$ $\alpha,\beta,x$ $\tilde{\alpha},\tilde{\beta},\tilde{x}$

Advertência sobre a não linearidade

$y,x$ $y$ $\tilde{x}$ $x$ $x$ $\tilde{x}$ $x$ $y$ $\tilde{x}$ $x$ $\mathbb{E}(x \mid \tilde{x})$ $\hat{\hat{y}}_x$ $x$

Código MATLAB para replicar o resultado do teste

Observe que isso também contém minhas próprias implementações para avaliar e OLS_solver, pois elas não foram fornecidas na pergunta.

rng(1)

OLS_solver = @(X,Y) [X ones(size(X))]'*[X ones(size(X))] \ ([X ones(size(X))]' * Y);
evaluate = @(b,x,y)  mean(([x ones(size(x))]*b - y).^2);

reg_mse_agg = [];
fixed_mse_agg = [];
reg_mse_orig_agg = [];
fixed_mse_orig_agg = [];

varMult = 1;
numTests = 60;
for dataNumber=1:8
    reg_mses = [];
    fixed_mses = [];
    reg_mses_orig = [];
    fixed_mses_orig = [];

    X = rand(1000,1);
    X(:,1) = X(:,1) * 10;
    X(:,1) = X(:,1) + 5;

    varX = var(X);
    y = 0.5 * X(:,1) -10;

    y = y + normrnd(0,1,size(y));    
    origX = X;
    X = X + normrnd(0,dataNumber * varMult ,size(X));
    train_size = floor(0.5 * length(y));
    for t=1:numTests,
        idx = randperm(length(y));
        train_idx = idx(1:train_size);
        test_idx = idx(train_size+1:end);
        Xtrain = X(train_idx,:);
        ytrain = y(train_idx);        
        Xtest = X(test_idx,:);
        origXtest = origX(test_idx,:);
        ytest = y(test_idx);

        b = OLS_solver(Xtrain, ytrain);
        %first arg of evaluate returns MSE, working correctly.
        reg_mse = evaluate( b,Xtest,ytest);
        reg_mses = [reg_mses ; reg_mse];

        varInd = var(Xtrain);
        varNoise = varInd - varX;
        bFixed = [0 0]';
        bFixed(1) = b(1) * varInd / varX;
        bFixed(2) = mean(ytrain - bFixed(1)*Xtrain);
        fixed_mse = evaluate( bFixed,Xtest,ytest);
        fixed_mses = [fixed_mses ; fixed_mse];

        reg_mse_orig = evaluate(b, origXtest, ytest);
        reg_mses_orig = [reg_mses; reg_mses_orig];

        fixed_mse_orig = evaluate(bFixed, origXtest, ytest);
        fixed_mses_orig = [fixed_mses_orig; fixed_mse_orig];

    end
    reg_mse_agg = [reg_mse_agg , reg_mses];
    fixed_mse_agg = [fixed_mse_agg , fixed_mses];

    reg_mse_orig_agg = [reg_mse_orig_agg , reg_mses_orig];
    fixed_mse_orig_agg = [fixed_mse_orig_agg , fixed_mses_orig]; 
end

disp('Reg parameters, noisy predictor')
disp(mean(reg_mse_agg))
disp('Fixed parameters, noisy predictor')
disp(mean(fixed_mse_agg))
disp('Reg parameters, true predictor')
disp(mean(reg_mse_orig_agg))
disp('Fixed parameters, true predictor')
disp(mean(fixed_mse_orig_agg))

Juho Kokkala
fonte