independent_metropolis_hastings_core.m 10.2 KB
Newer Older
1
function myoutput = independent_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
2
3
4
5
% PARALLEL CONTEXT
% The most computationally intensive portion of code in
% independent_metropolis_hastings (the 'for xxx = fblck:nblck' cycle).
% See the comment in random_walk_metropolis_hastings_core.m funtion.
6
7
%
% INPUTS
8
%   See See the comment in random_walk_metropolis_hastings_core.m funtion.
9

10
11
12
% OUTPUTS
%   See See the comment in random_walk_metropolis_hastings_core.m funtion.
%
13
14
% ALGORITHM
%   Portion of Independing Metropolis-Hastings.
15
16
17
18
%
% SPECIAL REQUIREMENTS.
%   None.
%
19
% Copyright (C) 2006-2011 Dynare Team
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%
% 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


41
global bayestopt_ estim_params_ options_  M_ oo_ objective_function_penalty_base
42

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
% 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;
60
InitSizeArray=myinputs.InitSizeArray;
61
62
63
64
65
record=myinputs.record;
varargin=myinputs.varargin;

if whoiam
    Parallel=myinputs.Parallel;
66
    % initialize persistent variables in priordens()
Marco Ratto's avatar
Marco Ratto committed
67
    priordens(xparam1',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
68
              bayestopt_.p3,bayestopt_.p4,1);
69
end
70

71
% (re)Set the penalty.
72
objective_function_penalty_base = Inf;
73

74
MhDirectoryName = CheckPath('metropolis',M_.dname);
75
76
77
78
79
80
81
82
83
84
85
86
87

OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
    n = npar;
    ProposalDensity = 'multivariate_normal_pdf';
elseif strcmpi(ProposalFun,'rand_multivariate_student')
    n = options_.student_degrees_of_freedom;
    ProposalDensity = 'multivariate_student_pdf';
end
% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
%%%%
%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
%%%%
88

89
90
91
92
93
94
95
96
97
98
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
99
100
101
102
103

jloop=0;

for b = fblck:nblck,
    jloop=jloop+1;
104
105
106
107
    try
        % this will not work if the master uses a random generator not
        % available in the slave (different Matlab version or
        % Matlab/Octave cluster). Therefor the trap.
108
109
110
111
112
113
        
        % this set the random generator type (the seed is useless but
        % needed by the function)
        set_dynare_seed(options_.DynareRandomStreams.algo,...
                        options_.DynareRandomStreams.seed);
        % this set the state 
114
115
116
        set_dynare_random_generator_state(record.Seeds(b).Unifor, ...
                                          record.Seeds(b).Normal);
    catch
117
118
        % if the state set by master is incompatible with the slave, we
        % only reseed 
119
120
        set_dynare_seed(options_.DynareRandomStreams.seed+b);
    end
121
    if (options_.load_mh_file~=0)  && (fline(b)>1) && OpenOldFile(b)
122
        load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
123
              '_blck' int2str(b) '.mat'])
124
125
126
127
128
129
130
        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
131
132
133
    if exist('OCTAVE_VERSION') || options_.console_mode
        diary off
        disp(' ')
134
135
136
    elseif whoiam
        %       keyboard;
        waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...'];
137
        %       waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
138
139
140
        if options_.parallel(ThisMatlab).Local,
            waitbarTitle=['Local '];
        else
141
            waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
142
        end
143
        fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
144
145
146
    else,
        hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
        set(hh,'Name','Metropolis-Hastings');
147

148
149
150
151
152
153
    end
    isux = 0;
    jsux = 0;
    irun = fline(b);
    j = 1;
    while j <= nruns(b)
154
        par = feval(ProposalFun, xparam1, proposal_covariance, n);
155
        if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
156
            try
157
                logpost = - feval(TargetFun, par(:),varargin{:});
158
159
160
161
162
163
164
            catch,
                logpost = -inf;
            end
        else
            logpost = -inf;
        end
        r = logpost - ilogpo2(b) + ...
165
166
            log(feval(ProposalDensity, ix2(b,:), xparam1, proposal_covariance, n)) - ...
            log(feval(ProposalDensity, par, xparam1, proposal_covariance, n));
167
168
169
        if (logpost > -inf) && (log(rand) < r)
            x2(irun,:) = par;
            ix2(b,:) = par;
170
            logpo2(irun) = logpost;
171
172
173
            ilogpo2(b) = logpost;
            isux = isux + 1;
            jsux = jsux + 1;
174
        else
175
176
177
178
            x2(irun,:) = ix2(b,:);
            logpo2(irun) = ilogpo2(b);
        end
        prtfrc = j/nruns(b);
179
180
181
        if exist('OCTAVE_VERSION') || options_.console_mode
            if mod(j, 10) == 0
                if exist('OCTAVE_VERSION')
182
183
184
                    if (whoiam==0),
                        printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
                    end
185
186
187
                else
                    fprintf('   MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acception rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
                end
188
            end
189
            if mod(j,50)==0 && whoiam,
190
191
                %             keyboard;
                waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
192
                fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab))
193
194
            end
        else
195
            if mod(j, 3)==0 && ~whoiam
196
                waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
197
            elseif mod(j,50)==0 && whoiam,
198
199
                %             keyboard;
                waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
200
                fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab))
201
202
            end
        end
203

204
        if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
            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,:);
232
                record.LastLogPost(b) = logpo2(end);
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
            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;
248
    if exist('OCTAVE_VERSION') || options_.console_mode
249
        printf('\n');
250
        diary on
251
252
253
    elseif ~whoiam
        close(hh);
    end
254
    [record.Seeds(b).Unifor, record.Seeds(b).Normal] = get_dynare_random_generator_state();
255
256
257
258
259
260
261
262
263
264
    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;