Code covered by the BSD License  

Highlights from
Toolkit on Econometrics and Economics Teaching

image thumbnail

Toolkit on Econometrics and Economics Teaching

by

 

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

ADAPT_REJECT_SAMPLE.m
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: 
% Adaptive rejection sampling
% 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
        
        % add points
        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
    disp('Adaptive rejection sampling finished')
    disp(['Number of draws = ',num2str(ndraws)])
    disp(['Points forming envelope = ',num2str(npoints)])
    disp(['Acceptance rate when take the first draw = ',num2str(accept_rate(1))])
    disp(['Acceptance rate when take the last  draw = ',num2str(accept_rate(end))])
    % scatter(1:ndraws, accept_rate),title('Acceptance rate in sampling'), xlabel('Draws')
end

Contact us