1function Scale = calibrate_mh_scale_parameter(ObjectiveFunction, CovarianceMatrix, Parameters, MhBounds, options, varargin) 2 3% Tune the MH scale parameter so that the overall acceptance ratio is close to AcceptanceTarget. 4% 5% INPUTS 6% - ObjectiveFunction [fhandle] Function (posterior kernel). 7% - CovarianceMatrix [double] n*n matrix, covariance matrix of the jumping distribution. 8% - Parameters [double] n*1 vector, parameter values. 9% - MhBounds [double] n*2 matrix, bounds on the possible values for the parameters. 10% - options [structure] content of options_.tune_mh_jscale. 11% - varargin [cell] Additional arguments to be passed to ObjectiveFunction. 12% 13% OUTPUTS 14% - Scale [double] scalar, optimal scale parameter for teh jumping distribution. 15 16% Fire up the wait bar 17hh = dyn_waitbar(0,'Tuning of the scale parameter...'); 18set(hh,'Name','Tuning of the scale parameter.'); 19 20% Intilialize various counters. 21j = 1; jj = 1; isux = 0; jsux = 0; i = 0; 22 23% Evaluate the objective function. 24logpo0 = - feval(ObjectiveFunction, Parameters, varargin{:}); 25logpo1 = logpo0; 26 27% Get the dimension of the problem. 28n = length(Parameters); 29 30% Initialize the correction on the scale factor. 31correction = 1.0; 32 33% Set the initial value of the scale parameter 34Scale = options.guess; 35 36% Transposition of some arrays. 37MhBounds = MhBounds'; 38Parameters = Parameters'; 39 40% Compute the Cholesky of the covariance matrix, return an error if the 41% matrix is not positive definite. 42try 43 dd = chol(CovarianceMatrix); 44catch 45 error('The covariance matrix has to be a symetric positive definite matrix!') 46end 47 48% Set parameters related to the proposal distribution 49if options.rwmh.proposal_distribution=='rand_multivariate_normal' 50 nu = n; 51elseif options.rwmh.proposal_distribution=='rand_multivariate_student' 52 nu = options.rwmh.student_degrees_of_freedom; 53end 54 55% Random Walk Metropolis Hastings iterations... 56while j<=options.maxiter 57 % Obtain a proposal (jump) 58 proposal = feval(options.rwmh.proposal_distribution, Parameters, Scale*dd, nu); 59 % If out of boundaries set the posterior kernel equal to minus infinity 60 % so that the proposal will be rejected with probability one. 61 if all(proposal > MhBounds(1,:)) && all(proposal < MhBounds(2,:)) 62 logpo0 = -feval(ObjectiveFunction, proposal(:), varargin{:}); 63 else 64 logpo0 = -inf; 65 end 66 % Move if the proposal is enough likely... 67 if logpo0>-inf && log(rand)<logpo0-logpo1 68 Parameters = proposal; 69 logpo1 = logpo0; 70 isux = isux + 1; 71 jsux = jsux + 1; 72 end% ... otherwise I don't move. 73 prtfrc = j/options.maxiter; 74 % Update the waitbar 75 if ~mod(j, 10) 76 dyn_waitbar(prtfrc, hh, sprintf('Acceptance ratio [during last 500]: %f [%f]', isux/j, jsux/jj)); 77 end 78 % Adjust the value of the scale parameter. 79 if ~mod(j, options.stepsize) 80 r1 = jsux/jj; % Local acceptance ratio 81 r2 = isux/j; % Overall acceptance ratio 82 % Set correction for the scale factor 83 c1 = r1/options.target; 84 if abs(c1-1)>.05 85 correction = correction^options.rho*c1^(1-options.rho); 86 else 87 correction = c1; 88 end 89 % Apply correction 90 if c1>0 91 Scale = Scale*correction; 92 else 93 Scale = Scale/10; 94 end 95 % Update some counters. 96 jsux = 0; jj = 0; 97 if abs(r2-options.target)<options.c2 && abs(r1-options.target)<options.c1 98 i = i+1; 99 else 100 i = 0; 101 end 102 % Test convergence. 103 if i>options.c3 104 break 105 end 106 end 107 j = j+1; 108 jj = jj + 1; 109end 110 111dyn_waitbar_close(hh);