Gibbs程序，很全面

function result = coda(draws,vnames,info,fid) % PURPOSE: MCMC convergence diagnostics, modeled after Splus coda % ------------------------------------------------------ % USAGE: coda(draws,vnames,info,fid) % or: result = coda(draws) % where: draws = a matrix of MCMC draws (ndraws x nvars) % vnames = (optional) string vector of variable names (nvar x 1) % info = (optional) structure setting input values % info.q = Raftery quantile (default = 0.025) % info.r = Raftery level of precision (default = 0.01) % info.s = Raferty probability for r (default = 0.950) % info.p1 = 1st % of sample for Geweke chi-sqr test (default = 0.2) % info.p2 = 2nd % of sample for Geweke chi-sqr test (default = 0.5) % fid = file id for printing to a file, e.g. fid = fopen('cout','w'); % ------------------------------------------------------ % NOTES: you may supply only some of the info-structure arguments % the remaining ones will take on default values % ------------------------------------------------------ % RETURNS: output to command window if nargout = 0 % autocorrelation estimates % Rafterty-Lewis MCMC diagnostics % Geweke NSE, RNE estimates % Geweke chi-sqr prob on means from info.p1 vs info.p2 % a results structure if nargout = 1 % result.ndraw = # of draws % result.nvar = # of variables % result.p1 = p1 input argument (or default) % result.p2 = p2 input argument (or default) % result.q = Raftery q-value % result.r = Raftery r-value % result.s = Raftery s-value % result(i).pmean = posterior mean for variable i % result(i).pstd = posterior std deviation % result(i).nse = nse assuming no serial correlation for variable i % result(i).rne = rne assuming no serial correlation for variable i % result(i).nse1 = nse using 4% autocovariance tapered estimate % result(i).rne1 = rne using 4% autocovariance taper % result(i).nse2 = nse using 8% autocovariance taper % result(i).rne2 = rne using 8% autocovariance taper % result(i).nse3 = nse using 15% autocovariance taper % result(i).rne3 = rne using 15% autocovariance taper % result(i).nburn = number of draws required for burn-in % 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(i).auto1 = autocorrelation at lag 1 % result(i).auto5 = autocorrelation at lag 5 % result(i).auto10= autocorrelation at lag 10 % result(i).auto50= autocorrelation at lag 50 % % ------------------------------------------------------- % SEE ALSO: gmoment, apm, raftery % ------------------------------------------------------- % REFERENCES: Geweke (1992), `Evaluating the accuracy of sampling-based % approaches to the calculation of posterior moments', in J.O. Berger, % J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of % the Fourth Valencia International Meeting on Bayesian Statistics, % pp. 169-194, Oxford University Press % Also: `Using simulation methods for Bayesian econometric models: % Inference, development and communication', at: www.econ.umn.edu/~bacc % Best, N.G., M.K. Cowles, and S.K. Vines (1995) CODA: Manual % version 0.30. Biostatistics Unit, Cambridge U.K. http://www.mrc-bsu.cam.ac.uk % ----------------------------------------------------------------- % written by: % James P. LeSage, Dept of Economics % Texas State University-San Marcos % 601 University Drive % San Marcos, TX 78666 % jlesage@spatial-econometrics.com num_draws = length(draws); if num_draws < 500 error('coda: at least 500 draws are required'); end; if nargout == 1 % don't print, return a structure pflag = 1; else pflag = 0; % print to the command window end; if nargin == 4 if ~isstruct(info) error('coda: must supply options as a structure'); end; nflag = 0; [vsize junk] = size(vnames); % user may supply a blank vnames argument if vsize > 0 nflag = 1; % we have variable names end; fields = fieldnames(info); nf = length(fields); q = 0.025; r = 0.01; s = 0.95; p1 = 0.2; p2 = 0.5; for i=1:nf if strcmp(fields{i},'q') q = info.q; elseif strcmp(fields{i},'r') r = info.r; elseif strcmp(fields{i},'s') s = info.s; elseif strcmp(fields{i},'p1') p1 = info.p1; elseif strcmp(fields{i},'p2'); p2 = info.p2; end; end; elseif nargin == 3 if ~isstruct(info) error('coda: must supply options as a structure'); end; nflag = 0; [vsize junk] = size(vnames); % user may supply a blank vnames argument if vsize > 0 nflag = 1; % we have variable names end; fields = fieldnames(info); nf = length(fields); q = 0.025; r = 0.01; s = 0.95; p1 = 0.2; p2 = 0.5; fid = 1; for i=1:nf if strcmp(fields{i},'q') q = info.q; elseif strcmp(fields{i},'r') r = info.r; elseif strcmp(fields{i},'s') s = info.s; elseif strcmp(fields{i},'p1') p1 = info.p1; elseif strcmp(fields{i},'p2'); p2 = info.p2; end; end; elseif nargin == 2 nflag = 1; % we have variable names q = 0.025; r = 0.01; s = 0.95; % set default values p1 = 0.2; p2 = 0.5; fid = 1; elseif nargin == 1 nflag = 0; % no variable names q = 0.025; r = 0.01; s = 0.95; % set default values p1 = 0.2; p2 = 0.5; fid = 1; else error('Wrong # of arguments to coda'); end; result.q = q; result.r = r; result.s = s; [ndraw nvar] = size(draws); if nflag == 0 % no variable names make some up Vname = []; for i=1:nvar Vname{i} = str2mat(['variable ',num2str(i)]); end; elseif (nflag == 1) % the user supplied variable names Vname = []; [tst_n nsize] = size(vnames); if tst_n ~= nvar error('Wrong # of variable names in coda -- check vnames argument'); end; nmax = min(nsize,16); % truncate vnames to 16-characters for i=1:nvar Vname{i} = vnames(i,1:nmax); end; end; % end of nflag issue % =======> do SACF diagnostics nlag = 50; aout = zeros(50,nvar); for i=1:nvar; aout(:,i) = sacf(draws(:,i),nlag,1); end; % pull out sacf's at 1,5,10,50 aprt = zeros(nvar,4); aprt(:,1) = aout(1,:)'; aprt(:,2) = aout(5,:)'; aprt(:,3) = aout(10,:)'; aprt(:,4) = aout(50,:)'; % ========> do Raftery-Lewis diagnostics rafout = raftery(draws,q,r,s); rout = zeros(nvar,5); for i=1:nvar; rout(i,1) = rafout(i).kthin; rout(i,2) = rafout(i).nburn; rout(i,3) = rafout(i).n; rout(i,4) = rafout(i).nmin; rout(i,5) = rafout(i).irl; end; % =========> do Geweke diagnostics geweke = momentg(draws); % =========> split sample into 1st p1 percent and last p2 percent % and run Geweke chi-squared test result(1).p1 = p1; result(1).p2 = p2; nobs1 = round(p1*ndraw); nobs2 = round(p2*ndraw); draws1 = draws(1:nobs1,:); draws2 = trimr(draws,nobs2,0); res1 = momentg(draws1); res2 = momentg(draws2); resapm = apm(res1,res2); if pflag == 0 % print results to command window % =======> print SACF diagnostics fprintf(1,'MCMC CONVERGENCE diagnostics \n'); fprintf(1,'Based on sample size = %10d \n',ndraw); fprintf(1,'Autocorrelations within each parameter chain \n'); vstring = 'Variable'; lstring1 = 'Lag 1'; lstring2 = 'Lag 5'; lstring3 = 'Lag 10'; lstring4 = 'Lag 50'; cnames = strvcat(lstring1,l