Mercurial > hg > ltpda

comparison m-toolbox/classes/@pest/mcmcPlot.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: simple tool that plots mcmc pest objects
5 %
6 % CALL: mcmcPlot(pest_obj,pl)
7 %
8 % Parameters: - pest_obj: pest object
9 % - pl: plist
10 %
11 % example: - mcmcPlot(p,plist('plotmatrix',true,'burnin',5000,'pdfs',true,'chain',[1 2 3 4 5 6]))
12 %
13 %<a href="matlab:utils.helper.displayMethodInfo('pest', 'mcmcPlot')">ParametersDescription</a>
14 %
15 % Nikos Oct 2011
16 %
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18
19 function varargout = mcmcPlot(varargin)
20
21 %%% Check if this is a call for parameters
22 if utils.helper.isinfocall(varargin{:})
23 varargout{1} = getInfo(varargin{3});
24 return
25 end
26
27 import utils.const.*
28 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
29
30 % Collect input variable names
31 in_names = cell(size(varargin));
32 for ii = 1:nargin,in_names{ii} = inputname(ii);end
33
34 % Collect all AOs and plists
35 [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
36 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
37
38 % Decide on a deep copy or a modify
39 p = copy(pests, nargout);
40
41 % combine plists
42 pl = parse(pl, getDefaultPlist());
43 BurnIn = find(pl, 'burnin');
44 nbins = find(pl, 'nbins');
45 paramarray = find(pl, 'chain');
46 %colorm = find(pl, 'colormap');
47
48 if ~all(isa(pests, 'pest'))
49 error('### mcmcPlot must be only applied to pest objects.');
50 end
51
52 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53
54 outfigs = [];
55 N = numel(p);
56
57 if (BurnIn == 1 && ((find(pl, 'plotmatrix'))))
58 utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '...
59 'better and more accurate display the burn-in section of the chains should be discarded.']));
60 elseif (BurnIn == 1 && ((find(pl, 'pdfs'))))
61 utils.helper.msg(msg.IMPORTANT, sprintf(['The burn-in field is left empty or equal to one. For '...
62 'better and more accurate display the burn-in section of the chains should be discarded.']));
63 end
64
65 for numpest=1:N
66
67 % compute PDF
68 chain=p(numpest).chain(:,2:size(p(numpest).chain,2));
69 p(numpest).computePdf(plist('BurnIn',BurnIn,'nbins',nbins));
70
71 if isempty(paramarray)
72 % plot chain field (skip 1st column where the Loglikelihood is stored)
73 outfigs = [outfigs ; figure];
74 data = plot(chain);
75 else
76 ch = p(numpest).chain(:,paramarray);
77 outfigs = [outfigs ; figure];
78 data = plot(ch);
79 end
80
81 if (find(pl, 'plotmatrix'));
82 chn = p(numpest).chain(BurnIn:size(chain,1),2:size(p(numpest).chain,2));
83 outfigs = [outfigs ; figure];
84 plotmatrix(chn);
85 end
86
87 if (find(pl, 'results'));
88 chainn = chain(BurnIn:size(chain,1),:);
89 utils.helper.msg(msg.IMPORTANT, sprintf('Results:'));
90 for ii = 1:(size(chainn,2))
91 mu = mean(chainn(:,ii));
92 sigma = std(chainn(:,ii));
93 res = [mu sigma];
94 utils.helper.msg(msg.IMPORTANT, sprintf(' %d \t',res));
95 end
96 end
97
98 if (find(pl, 'pdfs'));
99 outfigs = [outfigs ; figure];
100
101 if ~(find(pl, 'plotmatrix'))
102 chn = p.chain(BurnIn:size(p.chain(:,:),1),:);
103 end
104
105 a=p(numpest).pdf;
106 a(:,1) = [];
107 a(:,1) = [];
108
109 for kk =1:size(chn,2)
110 subplot(2,4,kk)
111 x = linspace(min(a(:,2*kk-1)),max(a(:,2*kk-1)),10);
112 h = bar(a(:,2*kk-1),a(:,2*kk));
113 hold on;
114 y=normpdf(x,mean(chn(:,kk)),std(chn(:,kk)));
115 s=sum(y);
116 y=y/s;
117 plot(x,y,'r-','LineWidth',2);
118 hold off;
119
120 shading interp % Needed to graduate colors
121
122 ch = get(h,'Children');
123 fvd = get(ch,'Faces');
124 fvcd = get(ch,'FaceVertexCData');
125 n = 10;
126 [zs, izs] = sortrows(a(:,2*kk),1);
127 k = 128; % Number of colors in color table
128 colormap(summer(k)); % Expand the previous colormap
129 shading interp % Needed to graduate colors
130 for i = 1:n
131 color = floor(k*i/n); % Interpolate a color index
132 row = izs(i); % Look up actual row # in data
133 fvcd(fvd(row,1)) = 1; % Color base vertices 1st index
134 fvcd(fvd(row,4)) = 1;
135 fvcd(fvd(row,2)) = color; % Assign top vertices color
136 fvcd(fvd(row,3)) = color;
137 end
138 set(ch,'FaceVertexCData', fvcd); % Apply the vertex coloring
139 set(ch,'EdgeColor','k') % Give bars black borders
140 end
141
142
143
144 end
145
146
147 end
148
149
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151
152 if nargout == 0
153 out = outfigs;
154 else
155 error('### mcmcPlot cannot be used as a modifier!');
156 end
157
158 % Set outputs
159 if nargout > 0
160 varargout{1} = out;
161 end
162
163 end
164
165
166 %--------------------------------------------------------------------------
167 % Get Info Object
168 %--------------------------------------------------------------------------
169 function ii = getInfo(varargin)
170 if nargin == 1 && strcmpi(varargin{1}, 'None')
171 sets = {};
172 pl = [];
173 else
174 sets = {'Default'};
175 pl = getDefaultPlist;
176 end
177 % Build info object
178 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);
179 end
180
181
182 %--------------------------------------------------------------------------
183 % Get Default Plist
184 %--------------------------------------------------------------------------
185 function plout = getDefaultPlist()
186 persistent pl;
187 if exist('pl', 'var')==0 || isempty(pl)
188 pl = buildplist();
189 end
190 plout = pl;
191 end
192
193 function pl = buildplist()
194 pl = plist();
195
196 p = param({'chain',['Insert an array containing the parameters to plot. If left empty,'...
197 'then by default will plot the chains of every parameter. (note: The loglikelihood is stored '...
198 'in the first column)']}, paramValue.DOUBLE_VALUE([]));
199 pl.append(p);
200
201 p = param({'BurnIn',['Number of samples (of the chains) to be discarded for the computation of the PDFs of the parameters. Also used'...
202 'for producing the plotmatrix figure.']}, paramValue.DOUBLE_VALUE(1));
203 pl.append(p);
204
205 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));
206 pl.append(p);
207
208 p = param({'plotmatrix','Boolean to determine if a plotmatrix is desired'}, {1, {false,true}, paramValue.OPTIONAL});
209 pl.append(p);
210
211 p = param({'pdfs','Boolean to determine if a plot of the PDFs of each parameter is desired'}, {1, {false,true}, paramValue.OPTIONAL});
212 pl.append(p);
213
214 %p = param({'colormap','Choose a default matlab colormap for the parameter histogarms.'}, paramValue.DOUBLE_VALUE(summer));
215 %pl.append(p);
216
217 p = param({'results',['Set to "true" if a table of the results of the estimated parameters is desired.'...
218 'The results are printed on screen in 2 columns: the 1st contains the mean value'....
219 'and the second the sigma. Burn-in field is requiered.']}, {1, {false,true}, paramValue.OPTIONAL});
220 pl.append(p);
221
222 end
223