[Dart-dev] [6831] DART/trunk/DART_LAB/matlab: have the rank histogram filter call a local norm_inv() routine
nancy at ucar.edu
nancy at ucar.edu
Wed Feb 26 17:17:48 MST 2014
Revision: 6831
Author: nancy
Date: 2014-02-26 17:17:47 -0700 (Wed, 26 Feb 2014)
Log Message:
-----------
have the rank histogram filter call a local norm_inv() routine
instead of the stats toolbox norminv().
have the oned_model call a local kurt() routine
instead of the stats toolbox kurtosis().
have plot_gaussian call a local norm_pdf() routine
instead of the stats toolbox normpdf().
have the run_lorenz_96 code call different increment
routines based on the menu button for EAKF, EnKF, or RHF.
add a new norm_inv() routine which computes the
inverse CDF for a normal distribution, either for
a single scalar or for an array. based on the
fortran code found in assim_tools_mod.f90.
with this update, DART_LAB should be completely
independent of the statistics toolbox.
Modified Paths:
--------------
DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
DART/trunk/DART_LAB/matlab/oned_model.m
DART/trunk/DART_LAB/matlab/plot_gaussian.m
DART/trunk/DART_LAB/matlab/run_lorenz_96.m
Added Paths:
-----------
DART/trunk/DART_LAB/matlab/norm_inv.m
-------------- next part --------------
Added: DART/trunk/DART_LAB/matlab/norm_inv.m
===================================================================
--- DART/trunk/DART_LAB/matlab/norm_inv.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/norm_inv.m 2014-02-27 00:17:47 UTC (rev 6831)
@@ -0,0 +1,79 @@
+function x = norm_inv(p)
+% Computes inverse CDF for the normal distribution,
+% evaluated at either the scalar value or array P.
+%
+% based on the code in DART assim_tools_mod.f90, which is
+% in turn based on http://home.online.no/~pjacklam/notes/invnorm
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% constants where precision is vital
+a1 = -39.69683028665376;
+a2 = 220.9460984245205;
+a3 = -275.9285104469687;
+a4 = 138.357751867269;
+a5 = -30.66479806614716;
+a6 = 2.506628277459239;
+b1 = -54.4760987982241;
+b2 = 161.5858368580409;
+b3 = -155.6989798598866;
+b4 = 66.80131188771972;
+b5 = -13.28068155288572;
+c1 = -0.007784894002430293;
+c2 = -0.3223964580411365;
+c3 = -2.400758277161838;
+c4 = -2.549732539343734;
+c5 = 4.374664141464968;
+c6 = 2.938163982698783;
+d1 = 0.007784695709041462;
+d2 = 0.3224671290700398;
+d3 = 2.445134137142996;
+d4 = 3.754408661907416;
+
+% Split into an inner and two outer regions which have separate fits
+p_low = 0.02425;
+p_high = 1 - p_low;
+
+% if P is scalar, simplify
+if (numel(p) == 1)
+ if (p < p_low)
+ q = sqrt(-2.0 * log(p));
+ x = (((((c1*q + c2)*q + c3)*q + c4)*q + c5)*q + c6) / ...
+ ((((d1*q + d2)*q + d3)*q + d4)*q + 1.0);
+ elseif (p >= p_low) & (p <= p_high)
+ q = p - 0.5;
+ r = q*q;
+ x = (((((a1*r + a2)*r + a3)*r + a4)*r + a5)*r + a6)*q / ...
+ (((((b1*r + b2)*r + b3)*r + b4)*r + b5)*r + 1.0);
+
+ else % (p > p_high)
+ q = sqrt(-2.0 * log(1.0 - p));
+ x = -(((((c1*q + c2)*q + c3)*q + c4)*q + c5)*q + c6) / ...
+ ((((d1*q + d2)*q + d3)*q + d4)*q + 1.0);
+ end
+
+else % if P is an array, so do it element by element
+
+ inds = (p < p_low);
+ q(inds) = sqrt(-2.0 * log(p(inds)));
+ x(inds) = (((((c1*q(inds) + c2).*q(inds) + c3).*q(inds) + c4).*q(inds) + c5).*q(inds) + c6) ./ ...
+ ((((d1*q(inds) + d2).*q(inds) + d3).*q(inds) + d4).*q(inds) + 1.0);
+
+ inds = (p >= p_low) & (p <= p_high);
+ q(inds) = p(inds) - 0.5;
+ r(inds) = q(inds).*q(inds);
+ x(inds) = (((((a1*r(inds) + a2).*r(inds) + a3).*r(inds) + a4).*r(inds) + a5).*r(inds) + a6).*q(inds) ./ ...
+ (((((b1*r(inds) + b2).*r(inds) + b3).*r(inds) + b4).*r(inds) + b5).*r(inds) + 1.0);
+
+ inds = (p > p_high);
+ q(inds) = sqrt(-2.0 * log(1.0 - p(inds)));
+ x(inds) = -(((((c1*q(inds) + c2).*q(inds) + c3).*q(inds) + c4).*q(inds) + c5).*q(inds) + c6) ./ ...
+ ((((d1*q(inds) + d2).*q(inds) + d3).*q(inds) + d4).*q(inds) + 1.0);
+end
+
+end
+
Property changes on: DART/trunk/DART_LAB/matlab/norm_inv.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
===================================================================
--- DART/trunk/DART_LAB/matlab/obs_increment_rhf.m 2014-02-26 21:33:32 UTC (rev 6830)
+++ DART/trunk/DART_LAB/matlab/obs_increment_rhf.m 2014-02-27 00:17:47 UTC (rev 6831)
@@ -165,7 +165,7 @@
np = p / alpha;
% Find spot in standard normal
-x = norminv(np, 0, 1);
+x = norm_inv(np);
% Add in the mean and normalize by sd
x = mean + x * sd;
Modified: DART/trunk/DART_LAB/matlab/oned_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/oned_model.m 2014-02-26 21:33:32 UTC (rev 6830)
+++ DART/trunk/DART_LAB/matlab/oned_model.m 2014-02-27 00:17:47 UTC (rev 6831)
@@ -152,7 +152,7 @@
ylabel('Kurtosis', 'FontSize', 14);
% Compute initial kurtosis
-handles.kurtosis = kurtosis(handles.ens);
+handles.kurtosis = kurt(handles.ens);
% Fourth and fifth subplot are for rank histogram
handles.r4 = subplot(4, 2, 7);
@@ -411,7 +411,7 @@
% Update moments
handles.error = abs(mean(handles.ens));
handles.spread = mean(abs(handles.ens - mean(handles.ens)));
-handles.kurtosis = kurtosis(handles.ens);
+handles.kurtosis = kurt(handles.ens);
% Update handles structure
guidata(hObject, handles);
@@ -529,7 +529,7 @@
ylabel('Kurtosis', 'FontSize', 14);
% Compute initial kurtosis
-handles.kurtosis = kurtosis(handles.ens);
+handles.kurtosis = kurt(handles.ens);
% Reset focus to the menu gui window
[gcbo_h, gcbo_fig] = gcbo;
@@ -638,7 +638,7 @@
% Plot the segement for the prior kurtosis
subplot(handles.r3);
- prior_kurtosis= kurtosis(ens_new);
+ prior_kurtosis= kurt(ens_new);
plot([handles.time_step - 1 + 0.1, handles.time_step - 0.1], ...
[handles.kurtosis, prior_kurtosis], 'r');
handles.kurtosis= prior_kurtosis;
@@ -652,9 +652,9 @@
% Update the rank data
subplot(handles.r4);
ens_rank = get_ens_rank(ens_new, 0);
- fprintf([sprintf('\ntimestep %d bin edges are ',handles.time_step), ...
- num2str(sort([-Inf ens_new Inf]),'%10.4f'),'\n'])
- fprintf('timestep %d bin/"rank" is %d\n',handles.time_step, ens_rank)
+ %fprintf([sprintf('\ntimestep %d bin edges are ',handles.time_step), ...
+ % num2str(sort([-Inf ens_new Inf]),'%10.4f'),'\n'])
+ %fprintf('timestep %d bin/"rank" is %d\n',handles.time_step, ens_rank)
% Plot the latest rank entry as a different color
temp_rank(:, 1) = handles.prior_rank(1:handles.ens_size + 1);
@@ -839,7 +839,7 @@
handles.spread= post_spread;
% Plot the segment for the updated kurtosis
- post_kurtosis= kurtosis(new_ens);
+ post_kurtosis= kurt(new_ens);
subplot(handles.r3);
plot([handles.time_step - 0.1, handles.time_step + 0.1], ...
[handles.kurtosis, post_kurtosis], 'r');
Modified: DART/trunk/DART_LAB/matlab/plot_gaussian.m
===================================================================
--- DART/trunk/DART_LAB/matlab/plot_gaussian.m 2014-02-26 21:33:32 UTC (rev 6830)
+++ DART/trunk/DART_LAB/matlab/plot_gaussian.m 2014-02-27 00:17:47 UTC (rev 6831)
@@ -15,7 +15,7 @@
num_points = 1001;
interval = x_range / num_points;
x = x_min:interval:x_max;
-y = weight * normpdf(x, mean, sd);
+y = weight * norm_pdf(x, mean, sd);
plot_handle = plot(x, y);
Modified: DART/trunk/DART_LAB/matlab/run_lorenz_96.m
===================================================================
--- DART/trunk/DART_LAB/matlab/run_lorenz_96.m 2014-02-26 21:33:32 UTC (rev 6830)
+++ DART/trunk/DART_LAB/matlab/run_lorenz_96.m 2014-02-27 00:17:47 UTC (rev 6831)
@@ -394,8 +394,19 @@
for i = 1:MODEL_SIZE
obs_prior = temp_ens(i, :);
obs(i) = handles.true_state(time, i) + obs_sd * randn;
+
% Compute the increments for observed variable
- [obs_increments, err] = obs_increment_eakf(obs_prior, obs(i), obs_error_var);
+ switch filter_type
+ case 'EAKF'
+ [obs_increments, err] = ...
+ obs_increment_eakf(obs_prior, obs(i), obs_error_var);
+ case 'EnKF'
+ [obs_increments, err] = ...
+ obs_increment_enkf(obs_prior, obs(i), obs_error_var);
+ case 'RHF'
+ [obs_increments, err] = ...
+ obs_increment_rhf(obs_prior, obs(i), obs_error_var);
+ end
% Regress the increments onto each of the state variables
for j = 1:MODEL_SIZE
More information about the Dart-dev
mailing list