modwt_wvar_ci_verification_tcase.m

时间序列分析中很用的源码,书的原名为时间序列分析的小波方法.

function tc = modwt_wvar_ci_verification_tcase
% modwt_wvar_ci_verification_tcase -- munit test case to test modwt_wvar_ci.
%
%****f*  wmtsa.Tests.dwt/modwt_wvar_ci_verification_tcase
%
% NAME
%   modwt_wvar_ci_verification_tcase -- munit test case to test modwt_wvar.
%
% USAGE
%   run_tcase('modwt_wvar_ci_verification_tcase')
%
% INPUTS
%
%
% OUTPUTS
%   tc            = tcase structure for modwt_wvar_ci_verification testcase.
%
% SIDE EFFECTS
%
%
% DESCRIPTION
%
%
% SEE ALSO
%   modwt_wvar_ci
%

% AUTHOR
%   Charlie Cornish
%
% CREATION DATE
%   2004-06-24
%
% COPYRIGHT
%
%
% CREDITS
%
%
% REVISION
%   $Revision: 612 $
%
%***

%   $Id: modwt_wvar_ci_verification_tcase.m 612 2005-10-28 21:42:24Z ccornish $

  
tc = MU_tcase_new(mfilename);

tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_edof_ocean_d6);
tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_nile_haar_622_715);
tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_nile_haar_716_1284);
tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_edof_ocean_d6_eta1);
tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_edof_ocean_d6_eta3);
tc = MU_tcase_add_test(tc, @test_verify_modwt_wvar_ci_16pt_haar_lv1);

return

function test_verify_modwt_wvar_ci_edof_ocean_d6(mode)
  % Test Description:  
  %    Verify edof and MJ results for presented in Table 333
  %    of ocean shear using d6 filter.

  [X, x_att] = wmtsa_data('msp_4096');

  N = length(X);
  wtfname = 'd6';
  J0 = 9;
  boundary = 'reflection';

  [WJt, VJt, w_att] = modwt(X, wtfname, J0, boundary);
  WJt = WJt(1:N,:);
  
  ci_method = 'chi2eta3';
  estimator = 'unbiased';

  % Calculate wvar
  [wvar, tmpCI_wvar, tmpedof, MJ] = ...
      modwt_wvar(WJt, 'none', estimator, wtfname);

  clear tmp*;

  wtf_s = modwt_filter(wtfname);
  ht = wtf_s.h;
  gt = wtf_s.g;
  L = wtf_s.L;

  LJ = equivalent_filter_width(L, 1:J0);
  NJt = modwt_num_nonboundary_coef(wtfname, N, 1:J0);
  lb = LJ;
  ub = ones(1, J0) * N;
  
  % Compute ci, eodf, Qeta, AJ

  [CI_wvar_eta1, edof_eta1, Qeta_eta1] = ...
      modwt_wvar_ci(wvar, NJt, 'chi2eta1', WJt, lb, ub);
  
  [CI_wvar_eta3, edof_eta3, Qeta_eta3] = ...
      modwt_wvar_ci(wvar, NJt, 'chi2eta3', WJt, lb, ub);

  [exp_wvar] = wmtsa_data('msp_d6_wvar');
  
  fuzzy_tolerance = 1E-15;
%  fuzzy_diff(exp_wvar, wvar, fuzzy_tolerance, 'summary')
  MU_assert_fuzzy_diff(exp_wvar, wvar, fuzzy_tolerance, ...
                       'wvars are not equal');

  % Check EDOFs of table 333.
  [exp_edof_eta1, exp_edof_eta2, exp_edof_eta3, exp_MJ] = ...
      wmtsa_data('table333');

  MU_assert_fuzzy_diff(exp_MJ, MJ, 0, ...
                       'MJs are not equal');


  MU_assert_fuzzy_diff(exp_edof_eta1, round(edof_eta1), 0, ...
                       'edof_eta1s are not equal');



  MU_assert_fuzzy_diff(exp_edof_eta3, round(edof_eta3), 0, ...
                       'edof_eta3s are not equal');

  
return

function test_verify_modwt_wvar_ci_nile_haar_622_715(mode)
  % Test Description:  
  %    Verify calculated modwt wvar to expected values.

  [X, x_att] = wmtsa_data('nile');

  X = X(1:94);
  N = length(X);

  wtfname = 'haar';
  J0 = 4;
  boundary = 'periodic';

  [WJt, VJt, w_att] = modwt(X, wtfname, J0, boundary);
    
  ci_method = 'chi2eta3';
  estimator = 'unbiased';

  % Calculate wvar
  [wvar, tmpCI_wvar, tmpedof, MJ] = ...
      modwt_wvar(WJt, 'none', estimator, wtfname);

  clear tmp*;

  wtf_s = modwt_filter(wtfname);
  ht = wtf_s.h;
  gt = wtf_s.g;
  L = wtf_s.L;

  LJ = equivalent_filter_width(L, 1:J0);
  NJt = modwt_num_nonboundary_coef(wtfname, N, 1:J0);
  lb = LJ;
  ub = ones(1, J0) * N;
  
  % Compute ci, eodf, Qeta, AJ
  [CI_wvar, edof, Qeta, AJ] = ...
      modwt_wvar_ci(wvar, NJt, ci_method, WJt, lb, ub);

  clear tmp*;

  % Check wvar

  [exp_wvar, exp_CI_wvar, exp_edof, exp_Qeta] = ...
      wmtsa_data('nile_wvar_622_715');

  % Check wvar
  fuzzy_tolerance = 1E-15;

  MU_assert_fuzzy_diff(exp_wvar, wvar, fuzzy_tolerance, ...
                       'wvars are not equal', mode);

  fuzzy_tolerance = 0;
  MU_assert_fuzzy_diff(exp_edof, edof, fuzzy_tolerance, ...
                       'edofs are not equal', mode);

  format long
  if (strcmp(mode, 'details'))
    exp_CI_wvar
    CI_wvar
  end
  
  
  fuzzy_tolerance = 1E-3;
  MU_assert_fuzzy_diff(exp_CI_wvar, CI_wvar, fuzzy_tolerance, ...
                       'CI_wvars are not equal', mode);

  if (strcmp(mode, 'details'))
    exp_Qeta
    Qeta
  end
  
  
  fuzzy_tolerance = 1E-6;
  MU_assert_fuzzy_diff(exp_Qeta, Qeta, fuzzy_tolerance, ...
                       'Qeta are not equal', mode);

return

function test_verify_modwt_wvar_ci_nile_haar_716_1284(mode)
  % Test Description:  
  %    Verify calculated modwt wvar to expected values.

  [X, x_att] = wmtsa_data('nile');

  X = X(95:end);
  N = length(X);

  wtfname = 'haar';
  J0 = 4;
  boundary = 'periodic';

  [WJt, VJt, att] = modwt(X, wtfname, J0, boundary);
    
  ci_method = 'chi2eta3';
  estimator = 'unbiased';

  % Calculate wvar
  [wvar, tmpCI_wvar, tmpedof, MJ] = ...
      modwt_wvar(WJt, 'none', estimator, wtfname);

  clear tmp*;

  wtf_s = modwt_filter(wtfname);
  ht = wtf_s.h;
  gt = wtf_s.g;
  L = wtf_s.L;

  LJ = equivalent_filter_width(L, 1:J0);
  NJt = modwt_num_nonboundary_coef(wtfname, N, 1:J0);
  lb = LJ;
  ub = ones(1, J0) * N;
  
  % Compute ci, eodf, Qeta, AJ
  [CI_wvar, edof, Qeta, AJ] = ...
      modwt_wvar_ci(wvar, NJt, ci_method, WJt, lb, ub);

  

  % Check wvar

  [exp_wvar, exp_CI_wvar, exp_edof, exp_Qeta] = ...
      wmtsa_data('nile_wvar_716_1284');

  
  format long
  
  fuzzy_tolerance = 1E-15;

  MU_assert_fuzzy_diff(exp_wvar, wvar, fuzzy_tolerance, ...
                       'wvars are not equal', mode);

  % Check CIs, EDOFs, Qetas
  [CI_wvar, edof, Qeta] = ...
      modwt_wvar_ci(wvar, MJ, ci_method, estimator, wtfname, WJt);

  
  fuzzy_tolerance = 0;
  MU_assert_fuzzy_diff(exp_edof, edof, fuzzy_tolerance, ...
                       'edofs are not equal', mode);

  format long
  if (strcmp(mode, 'details'))
    exp_CI_wvar
    CI_wvar
  end
  
  fuzzy_tolerance = 1E-3;
  MU_assert_fuzzy_diff(exp_CI_wvar, CI_wvar, fuzzy_tolerance, ...
                       'CI_wvars are not equal', mode);

  if (strcmp(mode, 'details'))
    exp_Qeta
    Qeta
  end

  fuzzy_tolerance = 1E-5;
  MU_assert_fuzzy_diff(exp_Qeta, Qeta, fuzzy_tolerance, ...
                       'Qeta are not equal', mode);

return

function test_verify_modwt_wvar_ci_edof_ocean_d6_eta1(mode)
  % Test Description:  
  %    Verify edof, CI_wvar, Qeta, AJ, MJ results for presented in Figure 333
  %    of ocean shear using d6 filter and chi-square eta1 method.

  [X, x_att] = wmtsa_data('msp_4096');

  N = length(X);
  wtfname = 'd6';
  J0 = 9;
  boundary = 'reflection';

  [WJt, VJ0t, att] = modwt(X, wtfname, J0, boundary);
  WJt = WJt(1:N,:);

  ci_method = 'chi2eta1';
  estimator = 'unbiased';

  % Calculate wvar
  [wvar, tmpCI_wvar, tmpedof, MJ] = ...
      modwt_wvar(WJt, 'none', estimator, wtfname);

  clear tmp*;

  wtf_s = modwt_filter(wtfname);
  ht = wtf_s.h;
  gt = wtf_s.g;
  L = wtf_s.L;

  LJ = equivalent_filter_width(L, 1:J0);
  NJt = modwt_num_nonboundary_coef(wtfname, N, 1:J0);
  lb = LJ;
  ub = ones(1, J0) * N;
  
  % Compute ci, eodf, Qeta, AJ
  [CI_wvar, edof, Qeta, AJ] = ...
      modwt_wvar_ci(wvar, NJt, ci_method, WJt, lb, ub);

  [exp_wvar, exp_CI_wvar, exp_edof, exp_Qeta, exp_MJ, exp_AJ] = ...
      wmtsa_data('msp_d6_wvar_edof1');


  fuzzy_tolerance = 1E-15;

  % Compare wvar's
  MU_assert_fuzzy_diff(exp_wvar, wvar, fuzzy_tolerance, ...
                       'wvars are not equal', mode);

  format long;

  if (strcmp(mode, 'details'))
    exp_MJ
    MJ
  end

  % Compare MJ's
  fuzzy_tolerance = 0;
  MU_assert_fuzzy_diff(exp_MJ, MJ, fuzzy_tolerance, ...
                       'MJs are not equal', mode);


  
  format long e;
  wvar4 = wvar .* wvar;
  exp_AJ_calc_from_exp_edof = (exp_MJ .* (exp_wvar .* exp_wvar)) ./ exp_edof;

  % Compare AJs