Code covered by the BSD License

# Toolkit on Econometrics and Economics Teaching

### Hang Qian (view profile)

Many MATLAB routines related to econometrics, statistics and introductory economics teaching.

```function sample = ADAPT_REJECT_SAMPLE(MyPDF_log, domain, MyDerivative, points, ndraws, npoints_max, diag_flag)

% Purpose:
% Generate random numbers from user-supplied p.d.f. (in logarithm)
% -----------------------------------
% Algorithm:
% Pair-wise linear function to form envelope of user-supplied p.d.f.
% -----------------------------------
% Usage:
% MyPDF_log = function handle of user-supplied p.d.f. (in logarithm)
% domain = domain of the function, say [-inf,inf],  [0,inf],  [a,b]
% MyDerivative = function handle of the derivative of user-supplied p.d.f. (in logarithm)
% points = a vector of points in the domain to form pair-wise linear function (finite real numbers)
% ndraws = number of draws intended to be taken
% npoints_max = points to form pair-wise envelope will automatically increase, set an upper limit
% diag_flag = whether to display diagnostics of sampling (1 = YES, 0 = NO)
% -----------------------------------
% Returns:
% sample = random numbers from user-supplied p.d.f (a vector of 1*ndraws)
% -----------------------------------
% Notes:
% 1. User-supplied p.d.f. (in logarithm) must be globally concave.
% 2. User-supplied p.d.f. needs to support vector inputs.
% 3. In user-supplied p.d.f., inputs outside domain should return 0 (-Inf for logarithm)
% 4. Points to form pair-wise envelope are all finite.
%
% Written by Hang Qian, Iowa State University
% Contact me:  matlabist@gmail.com

if nargin < 2;      domain = [-inf, inf];       end;
if nargin < 3;      MyDerivative = @(x) (MyPDF_log(x+1e-6) - MyPDF_log(x-1e-6)) / 2e-6; end
if nargin < 4;      points = [-1, 0, 1];    end
if nargin < 5 ;     ndraws = 1;   end
if nargin < 6;      npoints_max = 30;   end
if nargin < 7;      diag_flag = 1;      end

%--------------
% Diagnostics
%--------------
if diag_flag ==1
points(isinf(points)) = [];
points = sort(points);
end
%---------------------------------

% Config envelope of ln[f(x)]
% u(x) = a(i) * x + b(i) ,  k(i) < x < k(i+1)
sample = zeros(1,ndraws);
npoints = length(points);
points_diff = diff(points);
points_val = MyPDF_log(points);
a = MyDerivative(points);
b = points_val - a .* points;
k = - diff(b) ./ diff(a);
kbig = [domain(1), k, domain(2)];

%--------------
% Diagnostics
%--------------
if diag_flag == 1
if any(sort(a,'descend') ~= a)
error('User-supplied p.d.f. (in logarithm) is not globally concave.')
end
if domain(1) == -inf  &&  a(1) < 0
error('When lower bound of domain is negative infinity, derivative should be positive.')
end
if domain(2) == inf  &&  a(end) > 0
error('When upper bound of domain is infinity, derivative should be negative.')
end

ind0 = find( abs(a) < 1e-8 );
if ~isempty(ind0);
disp('Do not use the largest point in the domain.')
points(ind0) = points(ind0) + 1e-4;
disp(['Auto correction: ',num2str(points)])
disp(' ')
points_diff = diff(points);
points_val = MyPDF_log(points);
a = MyDerivative(points);
b = points_val - a .* points;
k = - diff(b) ./ diff(a);
kbig = [domain(1), k, domain(2)];
end

accept_rate = zeros(1,ndraws);
end
%---------------------------------

% exp[u(x)] as the source density
% tip_val is the function value at connectors
% BigConstant is a robust constant
% cell_prob is the prob at each pairwised cell
% cell_cum_prob is the cdf

tip_val_log = [a,a(end)] .* kbig + [b,b(end)];
BigConstant = max(tip_val_log);
tip_val = exp(tip_val_log - BigConstant);
ub_val = tip_val(2:end) ./ a;                   % *exp(BigConstant)
lb_val = tip_val(1:end-1) ./ a;                 % *exp(BigConstant)
cell_prob = ub_val - lb_val;                    % *exp(BigConstant)
constant = 1 / sum(cell_prob);               % *exp(-BigConstant)
cell_prob = constant * cell_prob;           % OK
cell_cum_prob = cumsum(cell_prob);   % OK
shift_val = cell_cum_prob - constant * ub_val;  %OK

ntrial = 0;
r = 1;
while r <= ndraws

ntrial = ntrial + 1;

% Sample from the source with inverse cdf
U = rand;
ind = find( cell_cum_prob > U , 1);
slope_use = a(ind);
intercept_use = b(ind);
cand_val = log(  slope_use * (U - shift_val(ind)) / constant  ) + BigConstant;
cand = ( cand_val - intercept_use ) / slope_use;

W = rand;
log_W = log(W);
% ind1 = find( kbig > cand , 1) - 1;
ind1 = ind;
%slope_use = a(ind1);
%intercept_use = b(ind1);
%source_val = slope_use * cand + intercept_use;
source_val = cand_val;

if  (slope_use > 0 && ( tip_val_log(ind1)   - source_val) > log_W) ...
|| (slope_use < 0 && ( tip_val_log(ind1+1)  - source_val) > log_W)

%--------------
% Diagnostics
%--------------
if diag_flag == 1 ; accept_rate(r) = 1 / ntrial;  ntrial = 0;   end
%--------------

sample(r) = cand;
r = r + 1;
continue
end

ind2 = find( points > cand , 1);
if ind2 >= 2

if  (slope_use > 0 && (points_val(ind2-1) - source_val) > log_W) ...
|| (slope_use < 0 && (points_val(ind2) - source_val) > log_W)

%--------------
% Diagnostics
%--------------
if diag_flag == 1;  accept_rate(r) = 1 / ntrial;  ntrial = 0; end
%---------------

sample(r) = cand;
r = r + 1;
continue
end

%weight = (cand - points(ind2-1)) / (points(ind2) - points(ind2-1));
weight = (cand - points(ind2-1)) / points_diff(ind2-1);
target_val_lb =  (1 - weight) * points_val(ind2-1) + weight * points_val(ind2);
if (target_val_lb - source_val) > log_W

%--------------
% Diagnostics
%--------------
if diag_flag == 1;  accept_rate(r) = 1 / ntrial;  ntrial = 0; end
%--------------

sample(r) = cand;
r = r + 1;
continue
end
end

% accept with prob exp[ ln[f(x*)] C u(x*) ]
target_val = MyPDF_log(cand);
if (target_val - source_val) > log_W

%--------------
% Diagnostics
%--------------
if diag_flag == 1;   accept_rate(r) = 1 / ntrial;  ntrial = 0;  end
%--------------

sample(r) = cand;
r = r + 1;
end

% Expand pair-wise envelope
if r < ndraws && npoints < npoints_max

new_a = MyDerivative(cand);
new_b = target_val - new_a * cand;
a = [a, new_a];
b = [b, new_b];
points = [points,cand];
points_val = [points_val,target_val];

% rearrange
[points, index] = sort(points);
points_diff = diff(points);
points_val = points_val(index);
a = a(index);
b = b(index);
k = - diff(b) ./ diff(a);
kbig = [domain(1), k, domain(2)];

tip_val_log = [a,a(end)] .* kbig + [b,b(end)];
BigConstant = max(tip_val_log);
tip_val = exp(tip_val_log - BigConstant);
ub_val = tip_val(2:end) ./ a;
lb_val = tip_val(1:end-1) ./ a;
cell_prob = ub_val - lb_val;
constant = 1 / sum(cell_prob);
cell_prob = constant * cell_prob;
cell_cum_prob = cumsum(cell_prob);
shift_val = cell_cum_prob - constant * ub_val;

npoints = npoints + 1;
end

end

%--------------
% Diagnostics
%--------------
if diag_flag == 1