From b0f4f2af78334a986a4e2f6be1954dcba74348ca Mon Sep 17 00:00:00 2001 From: Marco Ratto <marco.ratto@jrc.ec.europa.eu> Date: Thu, 16 Apr 2015 18:03:07 +0200 Subject: [PATCH] Harmonize use of front-end mcf_analysis.m in GSA and identification routines. --- matlab/gsa/filt_mc_.m | 33 +++++++++++++++++++++++++++------ matlab/plot_identification.m | 23 ++++++++++++++++++++--- 2 files changed, 47 insertions(+), 9 deletions(-) diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index 6d562206d3..4724801352 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -242,18 +242,39 @@ end if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only if options_.opt_gsa.pprior anam='rmse_prior_post'; + atitle='RMSE prior: Log Posterior Kernel'; else anam='rmse_mc_post'; + atitle='RMSE MC: Log Posterior Kernel'; end - stab_map_1(x, ipost(1:nfilt), ipost(nfilt+1:end), anam, 1,[],OutDir); - stab_map_2(x(ipost(1:nfilt),:),alpha2,pvalue,anam, OutDir); + + options_mcf.pvalue_ks = alpha; + options_mcf.pvalue_corr = pvalue; + options_mcf.alpha2 = alpha2; + options_mcf.param_names = char(bayestopt_.name); + options_mcf.fname_ = fname_; + options_mcf.OutputDirectoryName = OutDir; + options_mcf.amcf_name = anam; + options_mcf.amcf_title = atitle; + options_mcf.title = atitle; + options_mcf.beha_title = 'better posterior kernel'; + options_mcf.nobeha_title = 'worse posterior kernel'; + mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, options_); + if options_.opt_gsa.pprior - anam='rmse_prior_lik'; + anam = 'rmse_prior_lik'; + atitle = 'RMSE prior: Log Likelihood Kernel'; else anam='rmse_mc_lik'; - end - stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir); - stab_map_2(x(ilik(1:nfilt),:),alpha2,pvalue,anam, OutDir); + atitle = 'RMSE MC: Log Likelihood Kernel'; + end + options_mcf.amcf_name = anam; + options_mcf.amcf_title = atitle; + options_mcf.title = atitle; + options_mcf.beha_title = 'better likelihood'; + options_mcf.nobeha_title = 'worse likelihood'; + mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, options_); + else if options_.opt_gsa.ppost, rmse_txt=rmse_pmean; diff --git a/matlab/plot_identification.m b/matlab/plot_identification.m index 4f9fcf98ac..549224ef0d 100644 --- a/matlab/plot_identification.m +++ b/matlab/plot_identification.m @@ -247,13 +247,30 @@ else hist(log10(idelre.cond)) title('log10 of Condition number in the LRE model') dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_ident_COND' ],options_); + options_mcf.pvalue_ks = 0.1; + options_mcf.pvalue_corr = 0.001; + options_mcf.alpha2 = 0; + options_mcf.param_names = name; + options_mcf.fname_ = M_.fname; + options_mcf.OutputDirectoryName = IdentifDirectoryName; + options_mcf.beha_title = 'LOW condition nbr'; + options_mcf.nobeha_title = 'HIGH condition nbr'; + options_mcf.amcf_name = 'MC_HighestCondNumberLRE'; + options_mcf.amcf_title = 'MC Highest Condition Number LRE Model'; + options_mcf.title = 'MC Highest Condition Number LRE Model'; ncut=floor(SampleSize/10*9); [dum,is]=sort(idelre.cond); - [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberLRE', 1, [], IdentifDirectoryName, 0.1,'MC Highest Condition Number LRE Model'); + mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_); + options_mcf.amcf_name = 'MC_HighestCondNumberModel'; + options_mcf.amcf_title = 'MC Highest Condition Number Model Solution'; + options_mcf.title = 'MC Highest Condition Number Model Solution'; [dum,is]=sort(idemodel.cond); - [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberModel', 1, [], IdentifDirectoryName, 0.1,'MC Highest Condition Number Model Solution'); + mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_); + options_mcf.amcf_name = 'MC_HighestCondNumberMoments'; + options_mcf.amcf_title = 'MC Highest Condition Number Model Moments'; + options_mcf.title = 'MC Highest Condition Number Model Moments'; [dum,is]=sort(idemoments.cond); - [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberMoments', 1, [], IdentifDirectoryName, 0.1,'MC Highest Condition Number Model Moments'); + mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_); % [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName); % for j=1:nparam, % % ibeh=find(idemoments.Mco(j,:)<0.9); -- GitLab