matbugs_25feb06.m

matlab中实现openbugs或winbugs的功能调用

  fname = fld{fi};
  tmp = getfield(S(1), fname);
  sz = size(tmp);
  psz = prod(sz);
  data = zeros(C, psz);
  for c=1:C
    tmp = getfield(S(c), fname);
    %data = cat(1, data, tmp);
    data(c,:) = tmp(:)';
  end
  if sz(2) > 1 % vector or matrix variable
    data = reshape(data, [C sz]);
  end
  A = setfield(A, fname, data);
end
  
 
%%%%%%%%%%%%

function [Rhat, m, s] = EPSR(samples)
% Note - not currently used.
%
% function [R, m, s] = EPSR(samples)
% "estimated potential scale reduction" statistics due to Gelman and Rubin.
% samples(i,j) for sample i, chain j
% 
% R = measure of scale reduction - value below 1.1 means converged:
%                                  see Gelman p297
% m = mean(samples)
% s = std(samples)

% This is the same as the netlab function convcalc(samples')

[n m] = size(samples);
meanPerChain = mean(samples,1); % each column of samples is a chain
meanOverall = mean(meanPerChain);

% Rhat only works if more than one chain is specified.
if m > 1
  % between sequence variace
  B = (n/(m-1))*sum( (meanPerChain-meanOverall).^2);

  % within sequence variance
  varPerChain = var(samples);
  W = (1/m)*sum(varPerChain);

  vhat = ((n-1)/n)*W + (1/n)*B;
  Rhat = sqrt(vhat/(W+eps));
else
  Rhat = nan;
end
   
m = meanOverall;
s = std(samples(:));

%%%%%%%%%

function stats = computeStats(A)

fld = fieldnames(A);
N = length(fld);
stats = struct('Rhat',[], 'mean', [], 'std', []);
for fi=1:length(fld)
  fname = fld{fi};
  samples = getfield(A, fname);
  sz = size(samples);
  clear R m s
  % samples(c, s, i,j,k)
  Nchains = sz(1);
  Nsamples = sz(2);

  st_mean_per_chain = mean(samples, 2);
  st_mean_overall   = mean(st_mean_per_chain, 1);
  
  % "estimated potential scale reduction" statistics due to Gelman and
  % Rubin.
  if Nchains > 1
    B = (Nsamples/Nchains-1) * ...
       sum((st_mean_per_chain - repmat(st_mean_overall, [Nchains,1])).^2);
    varPerChain = var(samples, 0, 2);
    W = (1/Nchains) * sum(varPerChain);
    vhat = ((Nsamples-1)/Nsamples) * W + (1/Nsamples) * B;
    Rhat = sqrt(vhat./(W+eps));
  else
    Rhat = nan;
  end

  % reshape and take standard deviation over all samples, all chains
  samp_shape = size(squeeze(st_mean_overall));
  % padarray is here http://www.mathworks.com/access/helpdesk/help/toolbox/images/padarray.html
  %reshape_target = padarray(samp_shape, [0 1], Nchains * Nsamples, 'pre');
  reshape_target = [Nchains * Nsamples, samp_shape]; % fix from Andrew Jackson  [e-mail]
  reshaped_samples = reshape(samples, reshape_target);
  st_std_overall = std(reshaped_samples);

  if ~isnan(Rhat)
    stats.Rhat = setfield(stats.Rhat, fname, squeeze(Rhat));
  end
  
  % special case - if mean is a 1-d array, make sure it's long
  squ_mean_overall = squeeze(st_mean_overall);
  st_mean_size = size(squ_mean_overall);
  if (length(st_mean_size) == 2) && (st_mean_size(2) == 1)
    stats.mean = setfield(stats.mean, fname, squ_mean_overall');
  else
    stats.mean = setfield(stats.mean, fname, squ_mean_overall);
  end 
  
  stats.std = setfield(stats.std, fname, squeeze(st_std_overall));
end

%%%%%%%%%%%%

function S=bugs2mat(file_ind,file_out,dir)
%BUGS2MAT  Read (Win)BUGS CODA output to matlab structure
%
% S=bugs2mat(file_ind,file_out,dir)
%  file_ind - index file (in ascii format)
%  file_out - output file (in ascii format)
%  dir      - directory where the files are found (optional)
%  S        - matlab structure, with CODA variables as fields
%
% The samples are stored in added 1'st dimension,
% so that 2 x 3 variable R with 1000 samples would be
% returned as S.R(1000,2,3)
%
% Note1: the data is returned in a structure that makes extraction
% of individual sample sequencies easy: the sequencies are
% directly Nx1 double vectors, as for example S.R(:,1,2).
% The computed statistics must, however, be squeezed,
% as mean(S.R,1) is a 1x2x2 matrix.
%
% Note2: in variable names "." is replaced with "_"

% To change the output structure, edit the 'eval' line in the m-file.
% For example, to return all samples as a cell, wich possibly varying
% number of samples for elements of a multidimensional variable,
% cange the 'eval' line to
%    eval(['S.' varname '={samples};']);
% Then the samples of R(2,1) would be returned as cell S.R(2,1)

% (c) [e-mail], 2000
% 2003-01-14 [e-mail] - Replace "." with "_" in variable names
% slightly modified by Maryam Mahdaviani, August 2005 (to suppress redundant output)

if nargin>2,
  file_ind=[dir '/' file_ind];
  file_out=[dir '/' file_out];
end

ind=readfile(file_ind);

data=load(file_out);

Nvars=size(ind,1);
S=[];
for k=1:Nvars
  [varname,indexstr]=strtok(ind(k,:));
  varname=strrep(varname,'.','_');
  indices=str2num(indexstr);
  if size(indices)~=[1 2]
     error(['Cannot read line: [' ind(k,:) ']']);
  end
  sdata = size(data);
  %indices
  samples=data(indices(1):indices(2),2);
  varname(varname=='[')='(';
  varname(varname==']')=')';
  leftparen=find(varname=='(');
  outstruct=varname;
  if ~isempty(leftparen)
     outstruct=sprintf('%s(:,%s',varname(1:leftparen-1),varname(leftparen+1:end));
  end
  eval(['S.' outstruct '=samples;']);
end

function T=readfile(filename)
f=fopen(filename,'r');
if f==-1, fclose(f); error(filename); end
i=1;
while 1
 clear line;
 line=fgetl(f);
 if ~isstr(line), break, end
 n=length(line);
 T(i,1:n)=line(1:n);
 i=i+1;
end
fclose(f);


% PROCESS_OPTIONS - Processes options passed to a Matlab function.
%                   This function provides a simple means of
%                   parsing attribute-value options.  Each option is
%                   named by a unique string and is given a default
%                   value.
%
% Usage:  [var1, var2, ..., varn[, unused]] = ...
%           process_optons(args, ...
%                           str1, def1, str2, def2, ..., strn, defn)
%
% Arguments:   
%            args            - a cell array of input arguments, such
%                              as that provided by VARARGIN.  Its contents
%                              should alternate between strings and
%                              values.
%            str1, ..., strn - Strings that are associated with a 
%                              particular variable
%            def1, ..., defn - Default values returned if no option
%                              is supplied
%
% Returns:
%            var1, ..., varn - values to be assigned to variables
%            unused          - an optional cell array of those 
%                              string-value pairs that were unused;
%                              if this is not supplied, then a
%                              warning will be issued for each
%                              option in args that lacked a match.
%
% Examples:
%
% Suppose we wish to define a Matlab function 'func' that has
% required parameters x and y, and optional arguments 'u' and 'v'.
% With the definition
%
%   function y = func(x, y, varargin)
%
%     [u, v] = process_options(varargin, 'u', 0, 'v', 1);
%
% calling func(0, 1, 'v', 2) will assign 0 to x, 1 to y, 0 to u, and 2
% to v.  The parameter names are insensitive to case; calling 
% func(0, 1, 'V', 2) has the same effect.  The function call
% 
%   func(0, 1, 'u', 5, 'z', 2);
%
% will result in u having the value 5 and v having value 1, but
% will issue a warning that the 'z' option has not been used.  On
% the other hand, if func is defined as
%
%   function y = func(x, y, varargin)
%
%     [u, v, unused_args] = process_options(varargin, 'u', 0, 'v', 1);
%
% then the call func(0, 1, 'u', 5, 'z', 2) will yield no warning,
% and unused_args will have the value {'z', 2}.  This behaviour is
% useful for functions with options that invoke other functions
% with options; all options can be passed to the outer function and
% its unprocessed arguments can be passed to the inner function.

% Copyright (C) 2002 Mark A. Paskin
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
% USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [varargout] = process_options(args, varargin)

% Check the number of input arguments
n = length(varargin);
if (mod(n, 2))
  error('Each option must be a string/value pair.');
end

% Check the number of supplied output arguments
if (nargout < (n / 2))
  error('Insufficient number of output arguments given');
elseif (nargout == (n / 2))
  warn = 1;
  nout = n / 2;
else
  warn = 0;
  nout = n / 2 + 1;
end

% Set outputs to be defaults
varargout = cell(1, nout);
for i=2:2:n
  varargout{i/2} = varargin{i};
end

% Now process all arguments
nunused = 0;
for i=1:2:length(args)
  found = 0;
  for j=1:2:n
    if strcmpi(args{i}, varargin{j})
      varargout{(j + 1)/2} = args{i + 1};
      found = 1;
      break;
    end
  end
  if (~found)
    if (warn)
      warning(sprintf('Option ''%s'' not used.', args{i}));
      args{i}
    else
      nunused = nunused + 1;
      unused{2 * nunused - 1} = args{i};
      unused{2 * nunused} = args{i + 1};
    end
  end
end

% Assign the unused arguments
if (~warn)
  if (nunused)
    varargout{nout} = unused;
  else
    varargout{nout} = cell(0);
  end
end