Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

Notes on RegEM

View test_regem.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
% 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.