Skip to content

@pitakakariki /test_regem.m secret
Created

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
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
Something went wrong with that request. Please try again.