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);