Skip to content

GitLab

dynare
dynare
Commit d38c4de4 authored by Michel Juillard's avatar Michel Juillard
Browse files
Options

estimation_dll: adding test for logposterior

parent 667a25ce
Hide whitespace changes
Inline Side-by-side
Showing with 1436 additions and 0 deletions
mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m 0 → 100644
View file @ d38c4de4
function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% PARALLEL CONTEXT
% This function contain the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for'
% cycle and are completely independent than suitable to be executed in parallel way.
%
% INPUTS
% o myimput [struc] The mandatory variables for local/remote
% parallel computing obtained from random_walk_metropolis_hastings.m
% function.
% o fblck and nblck [integer] The Metropolis-Hastings chains.
% o whoiam [integer] In concurrent programming a modality to refer to the differents thread running in parallel is needed.
% The integer whoaim is the integer that
% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
% cluster.
% o ThisMatlab [integer] Allows us to distinguish between the
% 'main' matlab, the slave matlab worker, local matlab, remote matlab,
% ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS
% o myoutput [struc]
% If executed without parallel is the original output of 'for b =
% fblck:nblck' otherwise a portion of it computed on a specific core or
% remote machine. In this case:
% record;
% irun;
% NewFile;
% OutputFileName
%
% ALGORITHM
% Portion of Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS.
% None.
% PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed
% in parallel. The code sutable to be executed in parallel on multi core or cluster machine,
% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in
% parallel and called name_function.m and name_function_core.m.
% In addition in the parallel package we have second set of functions used
% to manage the parallel computation.
%
% This function was the first function to be parallelized, later other
% functions have been parallelized using the same methodology.
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management funtions.
% Copyright (C) 2006-2008,2010 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 <http://www.gnu.org/licenses/>.
if nargin<4,
whoiam=0;
end
global bayestopt_ estim_params_ options_ M_ oo_
% reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
TargetFun=myinputs.TargetFun;
ProposalFun=myinputs.ProposalFun;
xparam1=myinputs.xparam1;
vv=myinputs.vv;
mh_bounds=myinputs.mh_bounds;
ix2=myinputs.ix2;
ilogpo2=myinputs.ilogpo2;
ModelName=myinputs.ModelName;
fline=myinputs.fline;
npar=myinputs.npar;
nruns=myinputs.nruns;
NewFile=myinputs.NewFile;
MAX_nruns=myinputs.MAX_nruns;
d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
varargin=myinputs.varargin;
% Necessary only for remote computing!
if whoiam
Parallel=myinputs.Parallel;
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
bayestopt_.p3,bayestopt_.p4,1);
end
% (re)Set the penalty
bayestopt_.penalty = Inf;
MhDirectoryName = CheckPath('metropolis');
options_.lik_algo = 1;
OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
n = npar;
elseif strcmpi(ProposalFun,'rand_multivariate_student')
n = options_.student_degrees_of_freedom;
end
% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
%%%%
%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
%%%%
if any(isnan(bayestopt_.jscale))
if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6.
load([ModelName '_optimal_mh_scale_parameter'])
proposal_covariance = d*Scale;
else
error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.')
end
else
proposal_covariance = d*diag(bayestopt_.jscale);
end
jloop=0;
for b = fblck:nblck,
jloop=jloop+1;
randn('state',record.Seeds(b).Normal);
rand('state',record.Seeds(b).Unifor);
if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
'_blck' int2str(b) '.mat'])
x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
OpenOldFile(b) = 0;
else
x2 = zeros(InitSizeArray(b),npar);
logpo2 = zeros(InitSizeArray(b),1);
end
if exist('OCTAVE_VERSION') || options_.console_mode
diary off
disp(' ')
elseif whoiam
% keyboard;
waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...'];
% waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
if options_.parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
else,
hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
set(hh,'Name','Metropolis-Hastings');
end
isux = 0;
jsux = 0;
irun = fline(b);
j = 1;
while j <= nruns(b)
par = feval(ProposalFun, ix2(b,:), proposal_covariance, n);
if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) )
try
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
% testing logposterior DLL
[junk,logpost1] = logposterior(par(:),varargin{2},mexext);
if abs(logpost+logpost1) > 1e-10;
disp ([logpost -logpost1])
end
else
logpost = -inf;
end
if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b))
x2(irun,:) = par;
ix2(b,:) = par;
logpo2(irun) = logpost;
ilogpo2(b) = logpost;
isux = isux + 1;
jsux = jsux + 1;
else
x2(irun,:) = ix2(b,:);
logpo2(irun) = ilogpo2(b);
end
prtfrc = j/nruns(b);
if exist('OCTAVE_VERSION') || options_.console_mode
if mod(j, 10) == 0
if exist('OCTAVE_VERSION')
printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
else
fprintf(' MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acceptance rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
end
end
if mod(j,50)==0 & whoiam
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab));
end
else
if mod(j, 3)==0 & ~whoiam
waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
elseif mod(j,50)==0 & whoiam,
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
end
end
if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
fprintf(fidlog,[' Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
fprintf(fidlog,[' Posterior mean........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']);
fprintf(fidlog,[' Minimum value.........:\n']);;
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']);
fprintf(fidlog,[' Maximum value.........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
fprintf(fidlog,' \n');
fclose(fidlog);
jsux = 0;
if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:);
record.LastLogLiK(b) = logpo2(end);
end
% size of next file in chain b
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
% initialization of next file if necessary
if InitSizeArray(b)
x2 = zeros(InitSizeArray(b),npar);
logpo2 = zeros(InitSizeArray(b),1);
NewFile(b) = NewFile(b) + 1;
irun = 0;
end
end
j=j+1;
irun = irun + 1;
end% End of the simulations for one mh-block.
record.AcceptationRates(b) = isux/j;
if exist('OCTAVE_VERSION') || options_.console_mode
if exist('OCTAVE_VERSION')
printf('\n');
else
fprintf('\n');
end
diary on;
elseif ~whoiam
close(hh);
end
record.Seeds(b).Normal = randn('state');
record.Seeds(b).Unifor = rand('state');
OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
end% End of the loop over the mh-blocks.
myoutput.record = record;
myoutput.irun = irun;
myoutput.NewFile = NewFile;
myoutput.OutputFileName = OutputFileName;
\ No newline at end of file
mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m 0 → 100644
View file @ d38c4de4
C =[
-7.4073
-6.1860
-6.5983
-5.6088
-5.0547
-4.4774
-3.8081
-3.8425
-2.4178
-1.9835
-1.0395
-0.1583
-0.0397
0.3505
-0.1879
-0.0067
0.0478
-1.2247
-1.4349
-0.7973
-0.0461
0.5844
1.1372
1.3801
1.8023
2.2972
2.0469
2.5435
2.8169
3.2007
2.6705
3.0518
3.2445
3.8443
3.8525
4.9494
4.2770
4.9532
5.1441
3.7124
3.9880
3.6926
2.6005
1.8679
1.9085
1.5563
1.2308
0.3264
-0.2208
-0.2483
-0.4082
-1.0315
-1.6030
-1.5499
-1.3777
-2.1675
-2.5138
-2.8820
-2.6958
-2.4719
-1.9854
-1.7954
-2.2362
-1.0595
-0.8808
-0.8548
-1.2839
-0.1363
0.2104
0.8810
0.3555
0.4766
1.3269
1.4506
1.4308
1.6263
1.9842
2.3948
2.8710
3.0177
2.9305
3.1739
3.7380
3.8285
3.3342
3.7447
3.7830
3.1039
2.8413
3.0338
0.3669
0.0847
0.0104
0.2115
-0.6649
-0.9625
-0.7330
-0.8664
-1.4441
-1.0179
-1.2729
-1.9539
-1.4427
-2.0371
-1.9764
-2.5654
-2.8570
-2.5842
-3.0427
-2.8312
-2.3320
-2.2768
-2.1816
-2.1043
-1.8969
-2.2388
-2.1679
-2.1172
];
E =[
0.6263
0.7368
0.7477
1.0150
0.6934
0.4135
0.3845
0.2380
0.2853
0.5999
0.8622
1.2116
1.4921
1.5816
1.7259
1.6276
1.2422
0.8084
0.4710
-0.3704
-0.6427
-0.5323
-0.5562
-0.3651
-0.4356
-0.7164
-0.5816
-0.4635
-0.8456
-0.9708
-0.7138
-0.7499
-0.6941
-0.6656
-0.2912
-0.1650
0.0774
0.2307
0.4484
0.4942
0.4653
0.2196
0.1736
-0.1595
-0.3918
-0.4611
-0.8493
-0.7384
-1.0604
-1.2166
-1.7187
-1.6932
-1.7830
-1.7035
-2.2079
-2.3769
-2.2511
-2.1093
-2.4638
-2.4027
-2.1313
-1.9199
-1.7941
-1.4661
-1.2269
-1.0392
-1.0725
-0.7156
-0.4778
-0.4233
-0.0409
0.1620
0.4280
0.5873
1.0323
1.3420
1.6902
2.0680
2.8219
3.2511
3.2930
3.5633
3.8992
3.6874
3.2849
3.1614
2.6221
2.5067
1.9223
1.1777
0.4483
-0.0661
-0.4424
-0.9000
-1.1478
-1.2047
-1.1412
-1.2383
-1.1048
-0.9716
-0.9287
-1.0057
-1.0827
-1.0200
-1.0072
-1.1740
-1.2809
-1.1086
-0.9866
-0.8947
-0.5875
-0.2329
0.1493
0.4906
0.8400
1.0720
1.2648
1.5431
];
I =[
2.6617
2.4325
1.9592
3.2530
2.9949
3.7918
4.7444
4.8289
5.5983
7.8923
9.4297
9.5010
10.0150
10.0413
9.6046
6.4766
5.9647
3.0114
0.5683
-2.1226
-2.1855
-0.8329
-1.5207
-1.3419
-1.7897
-0.1476
0.4675
-1.6516
-1.5419
-1.3050
-1.2451
-0.7815
-0.7796
-0.3612
-2.4072
1.1162
1.1383
3.4132
5.0356
2.8016
2.1734
0.9366
-0.7050
-1.5021
-2.9868
-6.0237
-6.2589
-6.9138
-8.2340
-9.2589
-9.2465
-9.6988
-9.7782
-10.5645
-10.7544
-13.1583
-12.2718
-12.0131
-13.5983
-12.3579
-10.9146
-11.1572
-12.4935
-9.4393
-8.5535
-7.3723
-10.0169
-6.6088
-5.2045
-4.1024
-2.8472
-1.3139
0.0477
1.5629
3.6947
4.0327
4.1320
7.1400
9.1036
8.5609
7.6576
8.8022
8.9611
10.0871
9.4797
9.3964