Instantly share code, notes, and snippets.

# pitakakariki/test_regem.mSecret Created Oct 19, 2011

What would you like to do?
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')