raftery.m

计量工具箱

function  result = raftery(runs,q,r,s)
% PURPOSE: MATLAB version of Gibbsit by Raftery and Lewis (1991)
%    calculates the # of draws needed in MCMC to estimate
%    the posterior cdf of the q-quantile
%    to within +/- r accuracy with probability s
% --------------------------------------------------------------------
% USAGE:  raftery(draws,q,r,s)
% where:  draws = draws from the sampler (= ndraws x nvar matrix)
%             q = quantile of the quantity of interest
%             r = level of precision desired
%             s = probability associated with r
% --------------------------------------------------------------------        
% RETURNS: a structure
%     result.meth  = 'raftery'
%     result(i).nburn = number of draws required for burn-in (variable i)
%     result(i).nprec = number of draws required to achieve r precision
%     result(i).kthin = skip parameter for 1st-order Markov chain
%     result(i).irl   = I-statistic from Raftery and Lewis (1992)
%     result(i).kind  = skip parameter sufficient to get independence chain
%     result(i).nmin  = # draws if the chain is white noise
%     result(i).n     = nburn + nprec
%     result.draws    = # of draws in draws input matrix
%     result.nvar     = # of variables in draws input matrix
%     result.q        = q-value (from input)
%     result.r        = r-value (from input)
%     result.s        = s-value (from input)
% ---------------------------------------------------------------------
% NOTES:   Example values of q, r, s:                                 
%     0.025, 0.005,  0.95 (for a long-tailed distribution)             
%     0.025, 0.0125, 0.95 (for a short-tailed distribution);           
%     0.5, 0.05, 0.95;  0.975, 0.005, 0.95;  etc.                                                                                            
%  - The result is quite sensitive to r, being proportional to the       
%  inverse of r^2.                                                                                                                          
%  - For epsilon, we have always used 0.001.  It seems that the result   
%  is fairly insensitive to a change of even an order of magnitude in  
%  epsilon.                                                                                                                                  
%  - One way to use the program is to run it for several specifications  
%  of r, s and epsilon and vectors q on the same data set.  When one   
%  is sure that the distribution is fairly short-tailed, such as when  
%  q=0.025, then r=0.0125 seems sufficient.  However, if one is not    
%  prepared to assume this, safety seems to require a smaller value of 
%  r, such as 0.005.                                                                                                                        
%  - The program takes as input the name of a file containing an initial 
%  run from a MCMC sampler.  If the MCMC iterates are independent,     
%  then the minimum number required to achieve the specified accuracy  
%  is about $\Phi^{-1} (\frac{1}{2}(1+s))^2 q(1-q)/r^2$ and this would 
%  be a reasonable number to run first.                                
%  When q=0.025, r=0.005 and s=0.95, this number is 3,748;             
%  when q=0.025, r=0.0125 and s=0.95, it is 600.  
% ---------------------------------------------------------------------
% SEE ALSO: coda(), prt(), raftery_d.m                            
% ---------------------------------------------------------------------
% REFERENCES: Raftery, A.E.  and Lewis, S.M.  (1992).  How many iterations 
% in the  Gibbs sampler?  In Bayesian Statistics, Vol.  4 (J.M.  Bernardo, 
% J.O.   Berger, A.P.  Dawid and A.F.M.  Smith, eds.).  Oxford, U.K.: Oxford 
% University Press, 763-773.   This paper is available via the World Wide 
% Web by linking to URL  
% http://www.stat.washington.edu/tech.reports/pub/tech.reports % and then 
% selecting the "How Many Iterations in the Gibbs Sampler" % link.                                                              
%  This paper is also available via regular ftp using the following commands:                                                           
%    ftp ftp.stat.washington.edu (or 128.95.17.34)                                                                                                                                                   *
%  Raftery, A.E. and Lewis, S.M. (1992).  One long run with diagnos-   
%  tics: Implementation strategies for Markov chain Monte Carlo.       
%  Statistical Science, Vol. 7, 493-497.                                                                                                    *
%  Raftery, A.E. and Lewis, S.M. (1995).  The number of iterations,    
%  convergence diagnostics and generic Metropolis algorithms.  In      
%  Practical Markov Chain Monte Carlo (W.R. Gilks, D.J. Spiegelhalter  
%  and S. Richardson, eds.). London, U.K.: Chapman and Hall.           
%  This paper is available via the World Wide Web by linking to URL    
%    http://www.stat.washington.edu/tech.reports/pub/tech.reports      
%  and then selecting the "The Number of Iterations, Convergence       
%  Diagnostics and Generic Metropolis Algorithms" link.                
%  This paper is also available via ftp.stat.washington.edu (or 128.95.17.34)                               
% ----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% [e-mail]

% NOTE: this is a translation of FORTRAN code
%       (which is why it looks so bad)
 
 [n nvar] = size(runs); 
    tmp  = zeros(n,1); 
    work = zeros(n,1);
 
result.meth = 'raftery';
result.draws = n;
result.nvar = nvar;
result.q = q;
result.r = r;
result.s = s;
    
    for nv = 1:nvar; % big loop over variables
  
 if (q > 0.0),
   cutpt = empquant(runs,q);
   work = (runs <= cutpt);
 else,
  q = 0.0;
  i1 = find(runs == 0); i2 = find(runs == 1);
  ct1 = size(i1);   ct2 = size(i2); 
 if (ct1+ct2 ~= n), error('raftery needs 0s and 1s in runs'); end;
 work = runs;
 q = sum(runs);
 q = q/n;
 end;      % end of if;
 
 kthin = 1; bic = 1.0; epss = 0.001;
  
 while(bic > 0);
 [tcnt tmp] = thin(work,n,kthin);
 [g2 bic] = mctest(tmp,tcnt);
 kthin = kthin+1; 
 end;      % end of while
 
 kthin = kthin-1; 
 [alpha beta]=mcest(tmp,tcnt);
 kmind=kthin; 
 [g2 bic]=indtest(tmp,tcnt);
  
 while(bic > 0);
 [tcnt tmp] = thin(work,n,kmind); [g2 bic] = indtest(tmp,tcnt);
 kmind = kmind + 1; 
 end;      % end of while
  
 psum  = alpha + beta;
 tmp1  = log(psum*epss/max(alpha,beta)) / log(abs(1.0 - psum));
 nburn = fix((tmp1+1.0)*kthin);
 phi   = ppnd((s+1.0)/2);
 tmp2  = (2.0 - psum)*alpha*beta*(phi^2)/(psum^3 * r^2);
 nprec = fix(tmp2+1.0)*kthin; 
 nmin  = fix(((1.0-q)*q*phi^2/r^2)+1.0);
 irl   = (nburn + nprec)/nmin; 
 kind  = max(fix(irl+1.0),kmind);

result(nv).nburn = nburn;
result(nv).nprec = nprec;
result(nv).kthin = kthin;
result(nv).kind = kind;
result(nv).irl = irl;
result(nv).nmin = nmin;
result(nv).n = nburn+nprec;

end; % end of big loop over variables