Biodrowy estymator regresji osiągający lepsze wyniki niż obiektywny w modelu błędu w modelu zmiennych

Pracuję nad niektórymi danymi syntetycznymi dla modelu błędu w zmiennej dla niektórych badań. Obecnie mam pojedynczą zmienną niezależną i zakładam, że znam wariancję prawdziwej wartości zmiennej zależnej.

Dzięki tym informacjom mogę uzyskać obiektywny estymator dla współczynnika zmiennej zależnej.

Model:

$\tilde{x} = x + e_1$
$y = 0.5x -10 + e_2$
Gdzie: dla niektórych
$e_1\text{~}N(0,\sigma^2)$ $\sigma$
$e_2\text{~}N(0,1)$

Gdzie wartości są znane tylko dla każdej próbki, a także znane jest standardowe odchylenie rzeczywistej wartości dla próbki: . $y,\tilde{x}$ $x$ $\sigma_x$

Otrzymuję współczynnik stronniczości ( ) za pomocą OLS, a następnie dokonuję korekt za pomocą: $\hat{\beta}$

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

Widzę, że mój nowy, obiektywny estymator dla współczynnika jest znacznie lepszy (bliższy wartości rzeczywistej) z tym modelem, ale MSE pogarsza się niż używanie estymatora stronniczości.

Co się dzieje? Spodziewałem się, że wszechstronny estymator przyniesie lepsze wyniki niż stronniczy.

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

Wyniki:

tendencyjny MSE estymatora:

ans =

  Columns 1 through 7

    1.2171    1.6513    1.9989    2.3914    2.5766    2.6712    2.5997

  Column 8

    2.8346

MSE bezstronnego estymatora:

ans =

  Columns 1 through 7

    1.2308    2.0001    2.9555    4.9727    7.6757   11.3106   14.4283

  Column 8

   11.5653

Ponadto, drukowanie wartości bi bFixed- Widzę, że bFixedjest to rzeczywiście bliższe rzeczywistym wartościom 0.5,-10niż tendencyjny estymator (zgodnie z oczekiwaniami).

PS Wyniki bycia bezstronnym są gorsze od stronniczego estymatora są istotne statystycznie - test na to jest pominięty w kodzie, ponieważ jest to uproszczenie kodu „pełnej wersji”.

AKTUALIZACJA: Dodałem test, który sprawdza i , a stronniczy estymator jest rzeczywiście znacznie gorszy (większa wartość) niż obiektywny zgodnie z tą metryką, mimo że MSE tendencyjnego estymatora (w zestawie testowym) jest znacznie lepszy. Gdzie jest rzeczywistym współczynnikiem zmiennej zależnej, jest tendencyjnym estymatorem dla , a jest obiektywnym estymatorem dla . $\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$

To, jak sądzę, pokazuje, że przyczyną wyników NIE jest większa wariancja obiektywnego estymatora, ponieważ wciąż jest ona bliższa rzeczywistej wartości.

Źródło : Wykorzystanie notatek z wykładów Steve'a Pischke jako zasobu

— amit
źródło

Byłoby pomocne, gdybyś opublikował również wyniki, a nie tylko kod.

— Alecos Papadopoulos

@AlecosPapadopoulos Dodał go, nie dodał drukowania wszystkich wartości bi bFixed, ale wyjaśnił, co pokazują.

— amit

Podsumowanie: poprawione parametry służą do przewidywania jako funkcja prawdziwego predyktora . Jeśli w predykcji użyto , oryginalne parametry działają lepiej. $x$ $\tilde{x}$

Zauważ, że czają się dwa różne modele predykcji liniowej. Najpierw podane , drugie, podane , $y$ $x$

{\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}$

Nawet gdybyśmy mieli dostęp do prawdziwych parametrów, optymalne przewidywanie liniowe jako funkcje byłoby inne niż optymalne przewidywanie liniowe jako funkcja . Kod w pytaniu wykonuje następujące czynności $x$ $\tilde{x}$

Oszacuj parametry $\hat{\tilde{\beta}},\hat{\tilde{\alpha}}$
Oblicz szacunki $\hat{\beta},\hat{\alpha}$
Porównaj wydajność i $\hat{y}_1 = \hat{\beta}\,\tilde{x} + \hat{\alpha}$ $\hat{y}_2 = \hat{\tilde{\beta}}\,\tilde{x} + \hat{\tilde{\alpha}}$

Ponieważ w kroku 3 przewidujemy, że funkcja , a nie funkcja , lepsze jest użycie (szacunkowych) współczynników drugiego modelu. $\tilde{x}$ $x$

Rzeczywiście, gdybyśmy mieli dostęp do , i ale nie , moglibyśmy podstawić estymator liniowy w pierwszym modelu, Jeśli najpierw wykonamy transformację z postaci do a następnie wykonamy obliczenia w najnowszym równaniu, otrzymamy współczynniki $\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}$ . Tak więc, jeśli celem jest wykonanie prognozy liniowej, biorąc pod uwagę hałaśliwą wersję predyktora, powinniśmy po prostu dopasować model liniowy do hałaśliwych danych. Skorygowane współczynniki są interesujące, jeśli interesuje nas prawdziwe zjawisko z innych powodów niż przewidywanie.

α, β

$\alpha,\beta$

Testowanie

Edytowałem kod w OP, aby również oceniać MSE pod kątem prognoz, używając niehałasowej wersji prognozy (kod na końcu odpowiedzi). Wyniki są

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

Oznacza to, że gdy używasz zamiast , poprawione parametry rzeczywiście pokonują parametry nieskorygowane, zgodnie z oczekiwaniami. Co więcej, przewidywanie za pomocą ( ), to znaczy stałych parametrów i prawdziwego predyktora, jest lepsze niż ( ), że to parametry regu i hałaśliwy predyktor, ponieważ oczywiście hałas w pewnym stopniu szkodzi dokładności prognozowania. Pozostałe dwa przypadki odpowiadają użyciu parametrów niewłaściwego modelu, a zatem dają gorsze wyniki. $x$ $\tilde{x}$ $\alpha,\beta,x$ $\tilde{\alpha},\tilde{\beta},\tilde{x}$

Zastrzeżenie dotyczące nieliniowości

W rzeczywistości nawet jeśli związek między jest liniowy, związek między i może nie być. Zależy to od rozkładu . Na przykład w niniejszym kodzie jest pobierane z rozkładu jednolitego, więc bez względu na to, jak wysokie jest , znamy górną granicę dla a zatem przewidywane jako funkcja powinna nasycić. Możliwym rozwiązaniem w stylu bayesowskim byłoby ustalenie rozkładu prawdopodobieństwa dla a następnie podłączenie podczas wyprowadzania $y,x$ $y$ $\tilde{x}$ $x$ $x$ $\tilde{x}$ $x$ $y$ $\tilde{x}$ $x$ $\mathbb{E}(x \mid \tilde{x})$ $\hat{\hat{y}}_x$ - zamiast predykcji liniowej, której użyłem wcześniej. Jeśli jednak ktoś chce założyć rozkład prawdopodobieństwa dla , przypuszczam, że należy wybrać pełne rozwiązanie bayesowskie zamiast podejścia opartego przede wszystkim na korekcie oszacowań OLS. $x$

Kod MATLAB do replikacji wyniku testu

Zauważ, że zawiera to także moje własne implementacje do oceny i OLS_solver, ponieważ nie zostały podane w pytaniu.

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
źródło