Mercurial > hg > ltpda
view m-toolbox/classes/@pest/mcmcPlot.m @ 35:4be5f7a5316f database-connection-manager
Suppress output messages on preferences loading and saving
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DESCRIPTION: simple tool that plots mcmc pest objects % % CALL: mcmcPlot(pest_obj,pl) % % Parameters: - pest_obj: pest object % - pl: plist % % example: - mcmcPlot(p,plist('plotmatrix',true,'burnin',5000,'pdfs',true,'chain',[1 2 3 4 5 6])) % %<a href="matlab:utils.helper.displayMethodInfo('pest', 'mcmcPlot')">ParametersDescription</a> % % Nikos Oct 2011 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = mcmcPlot(varargin) %%% Check if this is a call for parameters if utils.helper.isinfocall(varargin{:}) varargout{1} = getInfo(varargin{3}); return end import utils.const.* utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename); % Collect input variable names in_names = cell(size(varargin)); for ii = 1:nargin,in_names{ii} = inputname(ii);end % Collect all AOs and plists [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names); pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); % Decide on a deep copy or a modify p = copy(pests, nargout); % combine plists pl = parse(pl, getDefaultPlist()); BurnIn = find(pl, 'burnin'); nbins = find(pl, 'nbins'); paramarray = find(pl, 'chain'); %colorm = find(pl, 'colormap'); if ~all(isa(pests, 'pest')) error('### mcmcPlot must be only applied to pest objects.'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% outfigs = []; N = numel(p); if (BurnIn == 1 && ((find(pl, 'plotmatrix')))) utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '... 'better and more accurate display the burn-in section of the chains should be discarded.'])); elseif (BurnIn == 1 && ((find(pl, 'pdfs')))) utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '... 'better and more accurate display the burn-in section of the chains should be discarded.'])); end for numpest=1:N % compute PDF chain=p(numpest).chain(:,2:size(p(numpest).chain,2)); p(numpest).computePdf(plist('BurnIn',BurnIn,'nbins',nbins)); if isempty(paramarray) % plot chain field (skip 1st column where the Loglikelihood is stored) outfigs = [outfigs ; figure]; data = plot(chain); else ch = p(numpest).chain(:,paramarray); outfigs = [outfigs ; figure]; data = plot(ch); end if (find(pl, 'plotmatrix')); chn = p(numpest).chain(BurnIn:size(chain,1),2:size(p(numpest).chain,2)); outfigs = [outfigs ; figure]; plotmatrix(chn); end if (find(pl, 'results')); chainn = chain(BurnIn:size(chain,1),:); utils.helper.msg(msg.IMPORTANT, sprintf('Results:')); for ii = 1:(size(chainn,2)) mu = mean(chainn(:,ii)); sigma = std(chainn(:,ii)); res = [mu sigma]; utils.helper.msg(msg.IMPORTANT, sprintf(' %d \t',res)); end end if (find(pl, 'pdfs')); outfigs = [outfigs ; figure]; if ~(find(pl, 'plotmatrix')) chn = p.chain(BurnIn:size(p.chain(:,:),1),:); end a=p(numpest).pdf; a(:,1) = []; a(:,1) = []; for kk =1:size(chn,2) subplot(2,4,kk) x = linspace(min(a(:,2*kk-1)),max(a(:,2*kk-1)),10); h = bar(a(:,2*kk-1),a(:,2*kk)); hold on; y=normpdf(x,mean(chn(:,kk)),std(chn(:,kk))); s=sum(y); y=y/s; plot(x,y,'r-','LineWidth',2); hold off; shading interp % Needed to graduate colors ch = get(h,'Children'); fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); n = 10; [zs, izs] = sortrows(a(:,2*kk),1); k = 128; % Number of colors in color table colormap(summer(k)); % Expand the previous colormap shading interp % Needed to graduate colors for i = 1:n color = floor(k*i/n); % Interpolate a color index row = izs(i); % Look up actual row # in data fvcd(fvd(row,1)) = 1; % Color base vertices 1st index fvcd(fvd(row,4)) = 1; fvcd(fvd(row,2)) = color; % Assign top vertices color fvcd(fvd(row,3)) = color; end set(ch,'FaceVertexCData', fvcd); % Apply the vertex coloring set(ch,'EdgeColor','k') % Give bars black borders end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargout == 0 out = outfigs; else error('### mcmcPlot cannot be used as a modifier!'); end % Set outputs if nargout > 0 varargout{1} = out; end end %-------------------------------------------------------------------------- % Get Info Object %-------------------------------------------------------------------------- function ii = getInfo(varargin) if nargin == 1 && strcmpi(varargin{1}, 'None') sets = {}; pl = []; else sets = {'Default'}; pl = getDefaultPlist; end % Build info object ii = minfo(mfilename, 'pest', 'ltpda', utils.const.categories.helper, '$Id: computePdf.m,v 1.2 2011/06/06 14:02:12 nikos Exp $', sets, pl); end %-------------------------------------------------------------------------- % Get Default Plist %-------------------------------------------------------------------------- function plout = getDefaultPlist() persistent pl; if exist('pl', 'var')==0 || isempty(pl) pl = buildplist(); end plout = pl; end function pl = buildplist() pl = plist(); p = param({'chain',['Insert an array containing the parameters to plot. If left empty,'... 'then by default will plot the chains of every parameter. (note: The loglikelihood is stored '... 'in the first column)']}, paramValue.DOUBLE_VALUE([])); pl.append(p); p = param({'BurnIn',['Number of samples (of the chains) to be discarded for the computation of the PDFs of the parameters. Also used'... 'for producing the plotmatrix figure.']}, paramValue.DOUBLE_VALUE(1)); pl.append(p); p = param({'nbins','Number of bins of the pdf histogram computed for every parameter (used again for the computation of the PDFs of the parameters)'}, paramValue.DOUBLE_VALUE(10)); pl.append(p); p = param({'plotmatrix','Boolean to determine if a plotmatrix is desired'}, {1, {false,true}, paramValue.OPTIONAL}); pl.append(p); p = param({'pdfs','Boolean to determine if a plot of the PDFs of each parameter is desired'}, {1, {false,true}, paramValue.OPTIONAL}); pl.append(p); %p = param({'colormap','Choose a default matlab colormap for the parameter histogarms.'}, paramValue.DOUBLE_VALUE(summer)); %pl.append(p); p = param({'results',['Set to "true" if a table of the results of the estimated parameters is desired.'... 'The results are printed on screen in 2 columns: the 1st contains the mean value'.... 'and the second the sigma. Burn-in field is requiered.']}, {1, {false,true}, paramValue.OPTIONAL}); pl.append(p); end