2006-07-05 16:06:45 +02:00
|
|
|
|
function metropolis(TargetFun,xparam1,vv,mh_bounds,varargin)
|
|
|
|
|
% Metropolis-Hastings algorithm.
|
|
|
|
|
%
|
|
|
|
|
% INPUTS
|
|
|
|
|
% o TargetFun [char] string specifying the name of the objective
|
|
|
|
|
% function (posterior kernel).
|
|
|
|
|
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
|
|
|
|
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
|
|
|
|
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
|
|
|
|
% o gend [integer] scalar specifying the number of observations ==> varargin{1}.
|
|
|
|
|
% o data [double] (T*n) matrix of data ==> varargin{2}.
|
|
|
|
|
%
|
|
|
|
|
% OUTPUTS
|
|
|
|
|
% None.
|
|
|
|
|
%
|
|
|
|
|
%
|
|
|
|
|
% ALGORITHM
|
|
|
|
|
% Metropolis-Hastings.
|
|
|
|
|
%
|
|
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
|
% None.
|
|
|
|
|
%
|
|
|
|
|
%
|
|
|
|
|
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
|
|
|
|
|
% Gnu Public License.
|
|
|
|
|
global M_ options_ bayestopt_
|
2005-02-18 20:54:39 +01:00
|
|
|
|
|
2006-01-16 13:58:18 +01:00
|
|
|
|
bayestopt_.penalty = 1e8;
|
|
|
|
|
|
2006-03-06 09:38:11 +01:00
|
|
|
|
MhDirectoryName = CheckPath('metropolis');
|
2005-02-18 20:54:39 +01:00
|
|
|
|
|
2005-09-10 12:04:22 +02:00
|
|
|
|
nblck = options_.mh_nblck;
|
|
|
|
|
nruns = ones(nblck,1)*options_.mh_replic;
|
|
|
|
|
npar = length(xparam1);
|
|
|
|
|
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
|
|
|
|
|
d = chol(vv);
|
|
|
|
|
options_.lik_algo = 1;
|
|
|
|
|
OpenOldFile = ones(nblck,1);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if options_.load_mh_file == 0
|
|
|
|
|
% Here we start a new metropolis-hastings, previous draws are not
|
|
|
|
|
% considered.
|
|
|
|
|
if nblck > 1
|
|
|
|
|
disp('MH: Multiple chains mode.')
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: One Chain mode.')
|
|
|
|
|
end
|
|
|
|
|
% Delete old mh files...
|
2006-03-06 09:38:11 +01:00
|
|
|
|
files = dir([ MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if length(files)
|
2006-03-06 09:38:11 +01:00
|
|
|
|
delete([ MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: Old _mh files succesfully erased!')
|
|
|
|
|
end
|
|
|
|
|
% Initial values...
|
|
|
|
|
if nblck > 1% Case 1: multiple chains
|
|
|
|
|
disp('MH: Searching for initial values...')
|
|
|
|
|
ix2 = zeros(nblck,npar);
|
|
|
|
|
ilogpo2 = zeros(nblck,1);
|
|
|
|
|
for j=1:nblck
|
|
|
|
|
validate = 0;
|
|
|
|
|
init_iter = 0;
|
|
|
|
|
trial = 1;
|
|
|
|
|
while validate == 0 & trial <= 10
|
|
|
|
|
candidate = options_.mh_init_scale*randn(1,npar)*d + transpose(xparam1);
|
|
|
|
|
if all(candidate' > mh_bounds(:,1)) & all(candidate' < mh_bounds(:,2))
|
|
|
|
|
ix2(j,:) = candidate;
|
2006-07-05 16:06:45 +02:00
|
|
|
|
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',varargin{:});
|
2005-09-10 12:04:22 +02:00
|
|
|
|
j = j+1;
|
|
|
|
|
validate = 1;
|
|
|
|
|
end
|
|
|
|
|
init_iter = init_iter + 1;
|
|
|
|
|
if init_iter > 100 & validate == 0
|
|
|
|
|
disp(['MH: I couldn''t get a valid initial value in 100 trials.'])
|
|
|
|
|
disp(['MH: You should Reduce mh_init_scale...'])
|
|
|
|
|
disp(sprintf('MH: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
|
|
|
|
|
options_.mh_init_scale = input('MH: Enter a new value... ');
|
|
|
|
|
trial = trial+1;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if trial > 10 & ~validate
|
|
|
|
|
disp(['MH: I''m unable to find a starting value for block ' int2str(j)])
|
|
|
|
|
return
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: Initial values found!')
|
|
|
|
|
disp(' ')
|
|
|
|
|
else% Case 2: one chain (we start from the posterior mode)
|
|
|
|
|
candidate = transpose(xparam1);
|
|
|
|
|
if all(candidate' > mh_bounds(:,1)) & all(candidate' < mh_bounds(:,2))
|
|
|
|
|
ix2 = candidate;
|
2006-07-05 16:06:45 +02:00
|
|
|
|
ilogpo2 = - feval(TargetFun,ix2',varargin{:});
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: Initialization at the posterior mode.')
|
2005-02-18 20:54:39 +01:00
|
|
|
|
disp(' ')
|
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: Initialization failed...')
|
|
|
|
|
disp('MH: The posterior mode lies outside the prior bounds.')
|
|
|
|
|
return
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
end
|
|
|
|
|
fblck = 1;
|
|
|
|
|
fline = ones(nblck,1);
|
|
|
|
|
NewFile = ones(nblck,1);
|
|
|
|
|
% Creation of the mh-history file:
|
2006-03-06 09:38:11 +01:00
|
|
|
|
file = dir([MhDirectoryName '/' M_.fname '_mh_history.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if length(files)
|
2006-03-06 09:38:11 +01:00
|
|
|
|
delete([ MhDirectoryName '/' M_.fname '_mh_history.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp('MH: Old mh_history file succesfully erased!')
|
|
|
|
|
end
|
|
|
|
|
AnticipatedNumberOfFiles = floor(nruns(1)/MAX_nruns);
|
|
|
|
|
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - AnticipatedNumberOfFiles*MAX_nruns;
|
|
|
|
|
record.Nblck = nblck;
|
|
|
|
|
record.MhDraws = zeros(1,3);
|
|
|
|
|
record.MhDraws(1,1) = nruns(1);
|
|
|
|
|
record.MhDraws(1,2) = AnticipatedNumberOfFiles+1;
|
|
|
|
|
record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
|
|
|
|
|
record.AcceptationRates = zeros(1,nblck);
|
|
|
|
|
record.Seeds.Normal = randn('state');
|
|
|
|
|
record.Seeds.Unifor = rand('state');
|
|
|
|
|
record.InitialParameters = ix2;
|
|
|
|
|
record.InitialLogLiK = ilogpo2;
|
|
|
|
|
record.LastParameters = zeros(nblck,npar);
|
|
|
|
|
record.LastLogLiK = zeros(nblck,1);
|
|
|
|
|
record.LastFileNumber = AnticipatedNumberOfFiles+1;
|
|
|
|
|
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
|
2006-03-06 11:00:32 +01:00
|
|
|
|
save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
|
2005-09-10 12:04:22 +02:00
|
|
|
|
elseif options_.load_mh_file == 1% Here we consider previous mh files (previous mh did not crash).
|
|
|
|
|
disp('MH: I''m loading past metropolis-hastings simulations...')
|
2006-03-06 09:38:11 +01:00
|
|
|
|
file = dir([ MhDirectoryName '/' M_.fname '_mh_history.mat' ]);
|
|
|
|
|
files = dir([ MhDirectoryName '/' M_.fname '_mh*.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if ~length(files)
|
|
|
|
|
disp('MH:: FAILURE! there is no MH file to load here!')
|
|
|
|
|
return
|
|
|
|
|
end
|
|
|
|
|
if ~length(file)
|
|
|
|
|
disp('MH:: FAILURE! there is no MH-history file!')
|
|
|
|
|
return
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2006-03-06 09:38:11 +01:00
|
|
|
|
load([ MhDirectoryName '/' M_.fname '_mh_history'])
|
2005-09-10 12:04:22 +02:00
|
|
|
|
end
|
|
|
|
|
past_number_of_blocks = record.Nblck;
|
|
|
|
|
if past_number_of_blocks ~= nblck
|
|
|
|
|
disp('MH:: The specified number of blocks doesn''t match with the previous number of blocks!')
|
|
|
|
|
disp(['MH:: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
|
|
|
|
|
disp(['MH:: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
|
|
|
|
|
nblck = past_number_of_blocks;
|
|
|
|
|
options_.mh_nblck = nblck;
|
|
|
|
|
end
|
|
|
|
|
% I get the last line of the last mh-file for initialization
|
|
|
|
|
% of the new metropolis-hastings simulations:
|
|
|
|
|
LastFileNumber = record.LastFileNumber;
|
|
|
|
|
LastLineNumber = record.LastLineNumber;
|
|
|
|
|
if LastLineNumber < MAX_nruns
|
|
|
|
|
NewFile = ones(nblck,1)*LastFileNumber;
|
|
|
|
|
else
|
|
|
|
|
NewFile = ones(nblck,1)*(LastFileNumber+1);
|
2006-07-05 16:06:45 +02:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
ilogpo2 = record.LastLogLiK;
|
|
|
|
|
ix2 = record.LastParameters;
|
|
|
|
|
fblck = 1;
|
|
|
|
|
fline = ones(nblck,1)*(LastLineNumber+1);
|
|
|
|
|
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
|
|
|
|
|
record.MhDraws = [record.MhDraws;zeros(1,3)];
|
|
|
|
|
NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
|
|
|
|
|
NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
|
|
|
|
|
AnticipatedNumberOfFiles = floor(NumberOfDrawsToBeSaved/MAX_nruns);
|
|
|
|
|
AnticipatedNumberOfLinesInTheLastFile = NumberOfDrawsToBeSaved - AnticipatedNumberOfFiles*MAX_nruns;
|
|
|
|
|
record.LastFileNumber = LastFileNumber + AnticipatedNumberOfFiles + 1;
|
|
|
|
|
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
|
|
|
|
|
record.MhDraws(end,1) = nruns(1);
|
|
|
|
|
record.MhDraws(end,2) = AnticipatedNumberOfFiles+1;
|
|
|
|
|
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
|
|
|
|
|
randn('state',record.Seeds.Normal);
|
|
|
|
|
rand('state',record.Seeds.Unifor);
|
2006-03-06 09:38:11 +01:00
|
|
|
|
save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
|
|
|
|
|
disp(' ')
|
|
|
|
|
elseif options_.load_mh_file == -1% The previous metropolis-hastings
|
|
|
|
|
% crashed before the end! I try to
|
|
|
|
|
% recover the saved draws...
|
|
|
|
|
disp('MH: Recover mode!')
|
|
|
|
|
disp(' ')
|
2006-03-06 09:38:11 +01:00
|
|
|
|
file = dir([MhDirectoryName '/' M_.fname '_mh_history.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if ~length(file)
|
|
|
|
|
disp('MH:: FAILURE! there is no MH-history file!')
|
|
|
|
|
return
|
|
|
|
|
else
|
2006-03-06 09:38:11 +01:00
|
|
|
|
load([ MhDirectoryName '/' M_.fname '_mh_history'])
|
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
nblck = record.Nblck;
|
|
|
|
|
options_.mh_nblck = nblck;
|
|
|
|
|
if size(record.MhDraws,1) == 1
|
|
|
|
|
OldMh = 0;% The crashed metropolis was the first session.
|
|
|
|
|
else
|
|
|
|
|
OldMh = 1;% The crashed metropolis wasn't the first session.
|
|
|
|
|
end
|
|
|
|
|
% Default initialization:
|
|
|
|
|
if OldMh
|
|
|
|
|
ilogpo2 = record.LastLogLiK;
|
|
|
|
|
ix2 = record.LastParameters;
|
|
|
|
|
else
|
|
|
|
|
ilogpo2 = record.InitialLogLiK;
|
|
|
|
|
ix2 = record.InitialParameters;
|
|
|
|
|
end
|
|
|
|
|
% Set "NewFile":
|
|
|
|
|
if OldMh
|
|
|
|
|
LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);
|
|
|
|
|
LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);
|
|
|
|
|
if LastLineNumberInThePreviousMh < MAX_nruns
|
|
|
|
|
NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
else
|
|
|
|
|
NewFile = ones(nblck,1);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
% Set fline (First line):
|
|
|
|
|
if OldMh
|
|
|
|
|
fline = ones(nblck,1)*(record.MhDraws(end-1,3)+1);
|
|
|
|
|
else
|
|
|
|
|
fline = ones(nblck,1);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
% Set fblck (First block):
|
|
|
|
|
fblck = 1;
|
|
|
|
|
% How many mh files should we have ?
|
|
|
|
|
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
|
|
|
|
|
ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
|
|
|
|
|
% I count the total number of saved mh files...
|
2006-03-06 09:38:11 +01:00
|
|
|
|
AllMhFiles = dir([MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
TotalNumberOfMhFiles = length(AllMhFiles);
|
|
|
|
|
% I count the number of saved mh files per block
|
|
|
|
|
NumberOfMhFilesPerBlock = zeros(nblck,1);
|
|
|
|
|
for i = 1:nblck
|
2006-03-06 09:38:11 +01:00
|
|
|
|
BlckMhFiles = dir([ MhDirectoryName '/' M_.fname '_mh*_blck' int2str(i) '.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
NumberOfMhFilesPerBlock(i) = length(BlckMhFiles);
|
|
|
|
|
end
|
|
|
|
|
tmp = NumberOfMhFilesPerBlock(1); b = 1;
|
|
|
|
|
% Is there a chain with less mh files than the first chain ?
|
|
|
|
|
CrashedBlck = 1;
|
|
|
|
|
while b <= nblck
|
|
|
|
|
if NumberOfMhFilesPerBlock(b) < ExpectedNumberOfMhFilesPerBlock
|
|
|
|
|
CrashedBlck = b;% YES!
|
|
|
|
|
break
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
b = b+1;
|
|
|
|
|
end
|
|
|
|
|
% The new metropolis-hastings should start from (fblck=CrashedBlck)
|
|
|
|
|
fblck = CrashedBlck;
|
|
|
|
|
% How many mh-files are saved in this block ?
|
|
|
|
|
NumberOfSavedMhFilesInTheCrashedBlck = ...
|
|
|
|
|
NumberOfMhFilesPerBlock(CrashedBlck);
|
|
|
|
|
% How many mh-files were saved in this block during the last session
|
|
|
|
|
% (if there was a complete session before the crash) ?
|
|
|
|
|
if OldMh
|
|
|
|
|
ante = sum(record.MhDraws(1:end-1,2),1);
|
2006-03-06 09:38:11 +01:00
|
|
|
|
load([MhDirectoryName '/' M_.fname '_mh' int2str(ante) '_blck' int2str(CrashedBlck) '.mat'],'logpo2');
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if length(logpo2) == MAX_nruns
|
|
|
|
|
IsTheLastFileOfThePreviousMhFull = 1;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
IsTheLastFileOfThePreviousMhFull = 0;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
else
|
|
|
|
|
ante = 0;% Because the crashed session is the first one
|
|
|
|
|
IsTheLastFileOfThePreviousMhFull = 1;
|
|
|
|
|
end
|
|
|
|
|
if ~IsTheLastFileOfThePreviousMhFull
|
|
|
|
|
MhFileExist = 1;
|
|
|
|
|
MhFileNumber = ante;
|
|
|
|
|
while MhFileExist
|
|
|
|
|
MhFileNumber = MhFileNumber + 1;
|
2006-03-06 09:38:11 +01:00
|
|
|
|
if ~exist([MhDirectoryName '/' M_.fname '_mh' int2str(MhFileNumber) '_blck' int2str(CrashedBlck) '.mat'])
|
2005-09-10 12:04:22 +02:00
|
|
|
|
MhFileExist = 0;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
% if MhFileNumber > ExpectedNumberOfMhFilesPerBlock % Peut <20>tre <20> d<>placer plus haut...
|
|
|
|
|
% disp('MH : You are using the recover mode but the previous session');
|
|
|
|
|
% disp(' of the metropolis-hastings didn''t crash!')
|
|
|
|
|
% disp('MH : I stop and you shoud modify your mod file.')
|
|
|
|
|
% return
|
|
|
|
|
% else
|
|
|
|
|
NumberOfCompletedMhFiles = (MhFileNumber-1)-ante;
|
|
|
|
|
% How many runs were saved ?
|
|
|
|
|
if OldMh
|
|
|
|
|
reste = MAX_nruns-record.MhDraws(end-1,3);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
reste = 0
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
NumberOfSavedDraws = MAX_nruns*(NumberOfCompletedMhFiles) + reste;
|
|
|
|
|
% Here is the number of draws we still need to complete the block:
|
|
|
|
|
nruns(CrashedBlck) = nruns(CrashedBlck)-NumberOfSavedDraws;
|
|
|
|
|
% I initialize with the last saved mh file of the inccomplete
|
|
|
|
|
% block:
|
2006-03-06 09:38:11 +01:00
|
|
|
|
load([MhDirectoryName '/' M_.fname '_mh' int2str(MhFileNumber-1) '_blck' int2str(CrashedBlck) '.mat']);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
ilogpo2(CrashedBlck) = logpo2(end);
|
|
|
|
|
ix2(CrashedBlck,:) = x2(end,:);
|
|
|
|
|
NewFile(CrashedBlck) = MhFileNumber;
|
|
|
|
|
fline(CrashedBlck,1) = 1;
|
|
|
|
|
else
|
|
|
|
|
% creuser ce qui se passe dans ce cas !
|
|
|
|
|
OpenOldFile(CrashedBlck) = 0;
|
|
|
|
|
disp('Ok')
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
end% of (if options_.load_mh_file == {0,1 or -1})
|
2005-02-18 20:54:39 +01:00
|
|
|
|
%%%%
|
2005-09-10 12:04:22 +02:00
|
|
|
|
%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
|
2005-02-18 20:54:39 +01:00
|
|
|
|
%%%%
|
2005-09-10 12:04:22 +02:00
|
|
|
|
InitSizeArray = min([MAX_nruns*ones(nblck) nruns],[],2);
|
|
|
|
|
for b = fblck:nblck
|
|
|
|
|
if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
|
2006-03-30 21:29:20 +02:00
|
|
|
|
load(['./' MhDirectoryName '/' M_.fname '_mh' int2str(NewFile(b)) ...
|
2006-04-03 16:36:24 +02:00
|
|
|
|
'_blck' int2str(b) '.mat'])
|
2005-09-10 12:04:22 +02:00
|
|
|
|
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);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(nblck) ')...']);
|
|
|
|
|
set(hh,'Name','Metropolis-Hastings')
|
|
|
|
|
isux = 0;
|
2006-07-05 16:06:45 +02:00
|
|
|
|
irun = fline(b);
|
2005-09-10 12:04:22 +02:00
|
|
|
|
j = 1;
|
|
|
|
|
while j <= nruns(b)
|
|
|
|
|
par = randn(1,npar)*d;
|
|
|
|
|
par = par.*bayestopt_.jscale' + ix2(b,:);
|
|
|
|
|
if all(par'>mh_bounds(:,1)) & all(par'<mh_bounds(:,2))
|
2006-07-05 16:06:45 +02:00
|
|
|
|
logpost = - feval(TargetFun,par',varargin{:});
|
2005-02-18 20:54:39 +01:00
|
|
|
|
else
|
2005-09-10 12:04:22 +02:00
|
|
|
|
logpost = -inf;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
if (logpost > -inf) & (log(rand) < logpost-ilogpo2(b))
|
|
|
|
|
x2(irun,:) = par;
|
|
|
|
|
ix2(b,:) = par;
|
|
|
|
|
logpo2(irun) = logpost;
|
|
|
|
|
ilogpo2(b) = logpost;
|
|
|
|
|
isux = isux + 1;
|
|
|
|
|
else
|
|
|
|
|
x2(irun,:) = ix2(b,:);
|
|
|
|
|
logpo2(irun) = ilogpo2(b);
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
prtfrc = j/nruns(b);
|
|
|
|
|
waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
|
|
|
|
|
if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
|
2006-03-06 09:38:11 +01:00
|
|
|
|
save([MhDirectoryName '/' M_.fname '_mh' int2str(NewFile(b)) '_blck' int2str(b)],'x2','logpo2');
|
2005-09-10 12:04:22 +02:00
|
|
|
|
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
|
|
|
|
|
if j == nruns(b) % I record the last draw...
|
|
|
|
|
record.LastParameters(b,:) = x2(end,:);
|
|
|
|
|
record.LastLogLiK(b) = logpo2(end);
|
|
|
|
|
end
|
|
|
|
|
if InitSizeArray(b)
|
|
|
|
|
x2 = zeros(InitSizeArray(b),npar);
|
|
|
|
|
logpo2 = zeros(InitSizeArray(b),1);
|
|
|
|
|
NewFile(b) = NewFile(b) + 1;
|
|
|
|
|
irun = 0;
|
|
|
|
|
else % InitSizeArray is equal to zero because we are at the end of an mc chain.
|
|
|
|
|
InitSizeArray(b) = min(nruns(b),MAX_nruns);
|
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2005-09-10 12:04:22 +02:00
|
|
|
|
j=j+1;
|
|
|
|
|
irun = irun + 1;
|
|
|
|
|
end% End of the simulations for one mh-block.
|
|
|
|
|
record.AcceptationRates(b) = isux/j;
|
|
|
|
|
close(hh);
|
|
|
|
|
end% End of the loop over the mh-blocks.
|
|
|
|
|
record.Seeds.Normal = randn('state');
|
|
|
|
|
record.Seeds.Unifor = rand('state');
|
2006-03-06 09:38:11 +01:00
|
|
|
|
save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
|
2005-09-10 12:04:22 +02:00
|
|
|
|
disp(['MH: Number of mh files : ' int2str(NewFile(1)) ' per block.'])
|
|
|
|
|
disp(['MH: Total number of generated files : ' int2str(NewFile(1)*nblck) '.'])
|
|
|
|
|
disp(['MH: Total number of iterations : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
|
2006-07-05 16:06:45 +02:00
|
|
|
|
disp(' ')
|