diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index 1b6e09ca2..aeca70b4e 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -268,6 +268,34 @@ if ~options_.opt_gsa.ppost [dum, ipost]=sort(-logpo2); [dum, ilik]=sort(-likelihood); end + +%% visual scatter analysis! +if options_.opt_gsa.ppost + tmp_title='R2 Posterior:'; + atitle='R2 Posterior:'; + asname='r2_post'; +else + if options_.opt_gsa.pprior + tmp_title='R2 Prior:'; + atitle='R2 Prior:'; + asname='r2_prior'; + else + tmp_title='R2 MC:'; + atitle='R2 MC:'; + asname='r2_mc'; + end +end +options_scatter.param_names = vvarvecm; +options_scatter.param_names_tex = vvarvecm_tex; +options_scatter.fname_ = fname_; +options_scatter.OutputDirectoryName = OutDir; +options_scatter.amcf_name = asname; +options_scatter.amcf_title = atitle; +options_scatter.title = tmp_title; +scatter_analysis(r2_MC, x,options_scatter, options_); +%% end of visual scatter analysis + + if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only if options_.opt_gsa.pprior anam='rmse_prior_post'; diff --git a/matlab/gsa/pick.m b/matlab/gsa/pick.m new file mode 100644 index 000000000..75c35927a --- /dev/null +++ b/matlab/gsa/pick.m @@ -0,0 +1,88 @@ +function pick + +% GLUEWIN, version 1.0, December 01. +% GLUEWIN is a MATLAB code designed for analysing the output +% of Monte Carlo runs when empirical observations of the model output are available +% and implements the GSA-GLUE methodology by Ratto et al. [1], based on a combination +% of GLUE (Generalised Likelihood Uncertainty Estimation) by K. Beven [2] and GSA +% Global Sensitivity Analysis) [3].'] +% The program has been developed by M. Ratto, European Commission, Joint Research Centre, +% Institute for the Protection and Security of The Citizen, Technological and Economic Risk Management, +% Applied Statistics, as a deliverable of the IMPACT project +% (EC Fifth Framework Programme, SCA Project, IST-1999-11313, DG-INFSO). +% +% The graphical layout of the code is inspired by the freeware GLUE package by K. Beven, +% vailable at the Lancaster University web site on the page [4]: +% http://www.es.lancs.ac.uk/hfdg/glue.html +% to which the GLUEWIN code introduces several extensions and additional options. +% Thanks are due to R. Girardi, A. Rossi, A. Saltelli, S. Tarantola and U. Callies for numerous +% comments and suggestions. +% For more information, please contact marco.ratto@jrc.it +% +% Disclaimer: This software has been developed at the Joint Research Centre of European Commission +% by officers in the course of their official duties. This software is not subject to copyright +% protection and is in the public domain. It is an experimental system. The Joint Research Centre +% of European Commission assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% +% [1] Ratto, M., Tarantola, S., A. Saltelli, Sensitivity analysis in model calibration: GSA-GLUE approach. +% 'Computer Physics Communications, 136, 2001, 212-224 +% [2] Beven K.J., Binley A., The Future of Distributed Models: Model Calibration and Uncertainty +% 'Prediction, Hydrological Processes, 6, 279-298, 1992 +% [3] Saltelli, A., K. Chan, M. Scott, Editors, (2000), Sensitivity analysis, John Wiley & Sons +% 'publishers, Probability and Statistics series. +% [4] Beven K., GLUE for Windows User manual, 1998. + + + +pmenu=findobj(gcf,'type','uicontextmenu','Tag','Run viewer'); +button1=findobj(gcf,'type','uimenu','Tag','save params'); +button2=findobj(gcf,'type','uimenu','Tag','eval params'); +%button=get(pmenu,'children'); +gg=gco; +ax0=gca; +set(gg,'buttondownfcn',[]); +c=get(gca,'currentpoint'); +x=c(1,1); +y=c(1,2); +X=get(gco,'xdata'); +Y=get(gco,'ydata'); +dx=get(gca,'xlim'); +dy=get(gca,'ylim'); +pos=get(gca,'position'); +scalex=dx(2)-dx(1); +scaley=dy(2)-dy(1); +if length(X)>1, + K = dsearchn([(Y./scaley)' (X./scalex)'],[y/scaley x/scalex]); +else + az=get(gca,'children'); + T =get(az(end),'ydata'); + [dum K]=max(T); +end + +KK=K; + +set(button1,'Label',['Save ',num2str(K)],'Callback',['scatter_callback(',num2str(KK),',''save'')']); +set(button2,'Label',['Eval ',num2str(K)],'Callback',['scatter_callback(',num2str(KK),',''eval'')']); +hh=findobj(gcf,'type','axes','Tag','scatter'); +for k=1:length(hh), + axes(hh(k)); + dum=get(gca,'children'); + dumx=get(dum(end),'xdata'); + dumy=get(dum(end),'ydata'); + xmid=min(dumx) + 0.5*(max(dumx)-min(dumx)); + hold on + plot(dumx(KK),dumy(KK),'or'); + if dumx(KK) < xmid, + text(dumx(KK),dumy(KK),[' ',num2str(K)], ... + 'FontWeight','Bold',... + 'Color','r'); + else + text(dumx(KK),dumy(KK),[num2str(K),' '], ... + 'HorizontalAlignment','right', ... + 'FontWeight','Bold',... + 'Color','r'); + end + hold off +end \ No newline at end of file diff --git a/matlab/gsa/scatter_analysis.m b/matlab/gsa/scatter_analysis.m new file mode 100644 index 000000000..15adfffd4 --- /dev/null +++ b/matlab/gsa/scatter_analysis.m @@ -0,0 +1,51 @@ +function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions) +% +% Written by Marco Ratto +% Joint Research Centre, The European Commission, +% marco.ratto@ec.europa.eu +% + +% Copyright (C) 2017 European Commission +% Copyright (C) 2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +param_names = options_scatter.param_names; + +if DynareOptions.TeX + if ~isfield(options_scatter,'param_names_tex') + param_names_tex = options_scatter.param_names; + else + param_names_tex = options_scatter.param_names_tex; + end +end +amcf_name = options_scatter.amcf_name; +amcf_title = options_scatter.amcf_title; +title = options_scatter.title; +fname_ = options_scatter.fname_; +xparam1=[]; +if isfield(options_scatter,'xparam1'), + xparam1=options_scatter.xparam1; +end +OutputDirectoryName = options_scatter.OutputDirectoryName; + +if ~DynareOptions.nograph, + skipline() + xx=[]; + if ~ isempty(xparam1), xx=xparam1; end + scatter_plots(lpmat, xdata, param_names, ... + '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, DynareOptions) +end diff --git a/matlab/gsa/scatter_callback.m b/matlab/gsa/scatter_callback.m new file mode 100644 index 000000000..8260cdf65 --- /dev/null +++ b/matlab/gsa/scatter_callback.m @@ -0,0 +1,42 @@ +function scatter_callback(K, type) +% +% Written by Marco Ratto +% Joint Research Centre, The European Commission, +% marco.ratto@ec.europa.eu +% + +% Copyright (C) 2017 European Commission +% Copyright (C) 2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +global oo_ M_ + +x=get(gcf,'userdata'); +r2=x{1}; +x=x{2}; + +xparam1=x(K,:)'; + +switch type + case 'save' + save(['my_params_' int2str(K)],'xparam1') + + case 'eval' + disp('Evaluating smoother ...') + oo_=evaluate_smoother(xparam1,M_.endo_names); + % [rmse, lnam, r2,vv] = plot_fit(obsname{:}); +end diff --git a/matlab/gsa/scatter_plots.m b/matlab/gsa/scatter_plots.m new file mode 100644 index 000000000..cbe9dfb32 --- /dev/null +++ b/matlab/gsa/scatter_plots.m @@ -0,0 +1,190 @@ +function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions) +% +% Written by Marco Ratto +% Joint Research Centre, The European Commission, +% marco.ratto@ec.europa.eu +% + +% Copyright (C) 2017 European Commission +% Copyright (C) 2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% PURPOSE: Pairwise scatter plots of the columns of x and y after +% Monte Carlo filtering +%--------------------------------------------------- +% USAGE: scatter_mcf(x,y,vnames,pltsym,diagon) +% or scatter_mcf(x,y) which relies on defaults +% where: +% x = an nxk matrix with columns containing behavioural sample +% y = an mxk matrix with columns containing non-behavioural sample +% vnames = a vector of variable names +% (default = numeric labels 1,2,3 etc.) +% pltsym = a plt symbol +% (default = '.' for npts > 100, 'o' for npts < 100 + + +[n,p] = size(X); +% X = X - ones(n,1)*min(Z); +% X = X ./ (ones(n,1)*max(Z)); + +nflag = 0; +if nargin >=3 + nflag = 1; +end; + +if nargin<4 || isempty(plotsymbol) + if n*p<100, plotsymbol = 'o'; + else plotsymbol = '.'; + end +end + +if nargin<5 + fnam=''; +end +if nargin<6, + dirname=''; + nograph=1; +else + nograph=0; +end +if nargin<7, + figtitle=fnam; +end +if nargin<8, + xparam1=[]; +end + +figtitle_tex=strrep(figtitle,'_','\_'); + +fig_nam_=[fnam]; +if ~nograph, + hh=dyn_figure(DynareOptions,'name',figtitle); + set(hh,'userdata',{X,xp}) +end + +bf = 0.1; +ffs = 0.05/(p-1); +ffl = (1-2*bf-0.05)/p; +if p>1, + fL = linspace(bf,1-bf+ffs,p+1); +else + fL = bf; +end +for i = 1:p + for j = 1:p + h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]); + if i==j + h1=cumplot(X(:,j)); + set(h,'Tag','cumplot') + % set(h1,'color',[0 0 1], 'linestyle','--','LineWidth',1.5) + set(h1,'color',[0 0 1],'LineWidth',1.5) + if ~isempty(xparam1) + hold on, plot(xparam1([j j]),[0 1],'k--') + end + if j

i + plot(X(:,i),X(:,j),[plotsymbol,'b']) + else + plot(X(:,i),X(:,j),[plotsymbol,'b']) + end + set(h,'Tag','scatter') + + %% + if ~isoctave + % Define a context menu; it is not attached to anything + hcmenu = uicontextmenu('Callback','pick','Tag','Run viewer'); + % Define callbacks for context menu + % items that change linestyle + hcb1 = 'scatter_callback'; + % hcb2 = ['set(gco,''LineStyle'','':'')']; + % hcb3 = ['set(gco,''LineStyle'',''-'')']; + % % Define the context menu items and install their callbacks + item1 = uimenu(hcmenu,'Label','save','Callback',hcb1,'Tag','save params'); + item2 = uimenu(hcmenu,'Label','eval','Callback',hcb1,'Tag','eval params'); + % item3 = uimenu(hcmenu,'Label','solid','Callback',hcb3); + % Locate line objects + hlines = findall(h,'Type','line'); + % Attach the context menu to each line + for line = 1:length(hlines) + set(hlines(line),'uicontextmenu',hcmenu) + end + end + %% + + if ~isempty(xparam1) + hold on, plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0]) + end + hold off; + % axis([-0.1 1.1 -0.1 1.1]) + if i