Create a gist now

Instantly share code, notes, and snippets.

Notes on RegEM
% 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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment