Notes on RegEM

test_regem.m

% s=randn('state');
randn('seed', 2)

T = cumsum(randn(500, 1)); % Random walk temperature (if not ok, what assumption that prevents this?)
R = randn(500, 6) * 10; % Random noise on proxies (if too much, what is the assumption for noise level?) 
Tm = T;

% 'RegEM assumes MAR, not MCAR' 
% 'All three cases are MAR'
% Trying cases 2 and 3:

obs = 1:200;
miss = 201:500;
% Tm(101:400,:)=NaN; % save rstate_middle s; % missing values at the middle (2), or
Tm(miss, :) = NaN; % save rstate_end s;  % missing values at the end (3). I didn't check the values itself

P = (repmat(T, 1, 6) + R) * diag(randn(6, 1)) * 10; % pseudoproxies with random scale

% Y= incomplete data matrix, which contains both instrumental
% data (surface temperature grid-box values
% arranged with rows representing the years and columns
% representing grid boxes) and proxy data 
Y=[Tm P];

opts = [];

[Xe, M, C, Xerr] = regem(Y, opts); % run the original regem without options

N = size(P, 1);

M1 = repmat(M(1), 1, N);
M2 = repmat(M(2:end), N, 1);
C11 = C(1, 1);
C12 = C(1, 2:end);
C22 = C(2:end, 2:end);

Tfit = M1 + C12 * (C22 \ (P - M2)');

Xerr2 = C11 - C12 * (C22 \ C12');
Xerr2 = Xerr2 * mean(Xerr(miss, 1) ./ Xerr2); % dof adjustment

mf_obs = mean(Tfit(obs));

lre = log(abs(Tfit(obs) - mf_obs) ./ abs(Tm(obs)' - mf_obs));
mre = exp(mean(lre));

Tfit_bc = mf_obs + (Tfit - mf_obs) / mre;

upper = mf_obs + ((Tfit - mf_obs) + 2 * Xerr2) / mre;
lower = mf_obs + ((Tfit - mf_obs) - 2 * Xerr2) / mre;

close all
plot(T,'r') % True
hold
plot(Xe(:,1), 'b') % Imputed
plot(Tfit, 'g') % Fitted
plot(upper, 'y') % CI Upper Bound
plot(lower, 'y') % CI Lower Bound

legend('True', 'Imputed', 'Fitted', 'CI Upper', 'CI Lower')