matbugs_25feb06.m.svn-base

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

function [samples, stats, structArray] = matbugs(dataStruct, bugsModel, varargin)
% MATBUGS a Matlab interface for WinBugs, similar to R2WinBUGS
% [samples, stats] = matbugs(dataStruct,  bugsModelFileName, ...)
%
% This generates a file called 'script.txt' in the directory where
% winbugs14.exe is kept (so you must have write access to this directory!),
% then calls winbugs with this script, then reads the resulting samples from
% files into matlab structs.
%
% INPUT:
% dataStruct contains values of observed variables.
% bugsModel is the name of the model file; MUST END IN .txt
% 
% Note: variables with names 'a.b' in the bugs model file
% should be called 'a_b' in the matlab data structure.
%
% Optional arguments passed as 'string', value pairs [default in brackets, case
% insensitive]:
%
% 'init' - init(i).v is a struct containing initial values for variable 'v'
%          for chain i.  Uninitialized variables are given random initial
%          values by WinBUGS.  It is highly recommended that you specify the
%          initial values of all root (founder) nodes.
% 'monitorParams' - cell array of field names (use 'a_b' instead of 'a.b')
%                   [defaults to *, which currently does nothing...]
% 'nChains'  - number of chains [3]
% 'nBurnin'  - num samples for burn-in per chain [1000]
% 'nSamples' - num samples to keep after burn-in [5000]
% 'thin'     - keep every n'th step [1]
% 'blocking' - do block updating for GLMs (using IRLS proposal)?
%              [1 - doesn't yet work]
% 'view'     - set to 1 if you want to view WinBUGS output (then close the WinBUGS
%                window to return control to matlab)
%              set to 0 to close WinBUGS automatically without pausing [default 0]
% 'Bugdir'   - location of winbugs executable ['C:/Program Files/WinBUGS14']
% 'workingDir' - directory to store temporary data/init/coda files [pwd/tmp]
%
% Note that total number of iterations  = nBurnin + nSamples * thin.
%
% OUTPUT
% S contains the samples; each field may have a different shape:
%  S.theta(c, s)       is the value of theta in sample s, chain c
%                      (scalar variable)
%  S.theta(c, s, i)    is the value of theta(i) in sample s, chain c
%                      (vector variable)
%  S.theta(c, s, i, j) is the value of theta(i,j) in sample s, chain c
%                      (matrix variable)
%
% stats contains various statistics, currently:
%    stats.mean, stats.std and stats.Rhat.
% Each field may have a different shape:
%    stats.mean.theta
%    stats.mean.theta(i)
%    state.mean.theta(i,j)
%
% Rhat is the "estimated potential scale reduction" statistic due to
%     Gelman and Rubin.
%
% Rhat values less than 1.1 mean the chain has probably converged for
%     this variable.
%
% Example
%
%     [S,stats] = matbugs(data, 'school_model.txt', 'nSamples', 10000, ..
%                          'monitorParams', {'mu_theta', 'theta'}, ...
%                          'init', initStruct, 'nchains', 3)
%
% Written by Maryam Mahdaviani ([e-mail])
% and Kevin Murphy ([e-mail]), August 2005
%
% Last updated 19 December 2005.


[initStructs, Bugdir, nChains, view, workingDir, nBurnin, nSamples, monitorParams, thin, blocking, junk] =  ...
   process_options(...
	varargin, ...
	'init', {}, ...
	'Bugdir', 'C:/Program Files/WinBUGS14', ...
	'nChains', 3, ...
	'view', 0, ...
	'workingDir', fullfile(pwd,'tmp'), ...
	'nBurnin', 1000, ...
	'nSamples', 5000, ...
	'monitorParams', {}, ...
	'thin', 1, ...
	'blocking', 1);

if length(initStructs) ~= nChains
  error(['init structure does not match number of chains ', ...
    sprintf('(%d)', nChains)]);
end

if ~exist(workingDir, 'dir')
  mkdir(pwd, 'tmp');
end

log_filename = fullfileKPM(workingDir, 'log.txt');
his_filename = fullfileKPM(workingDir, 'history.txt');

scriptFile = [Bugdir,'\','script.txt'];

%bugsModel = [pwd, '/', bugsModel];
%bugsModel = fullfile(pwd, bugsModel);
bugsModel = strrep(bugsModel, '\', '/'); % winBUGS wants a/b not a\b

codaFile = fullfileKPM(workingDir, 'coda');

fid = fopen(scriptFile,'w');
if (fid == -1)
   error(['Cannot open ', scriptFile]);
end
%display(option)
fprintf(fid, 'display(''log'') \n');

%check(model file)
fprintf(fid, 'check(''%s'')\n',bugsModel);

if ~isempty(dataStruct)
  dataFile = fullfileKPM(workingDir, 'data.txt');
  dataGen(dataStruct, dataFile);
  fprintf(fid, 'data(''%s'')\n', dataFile);
end
%fprintf(fid, '.txt'')\n');

fprintf(fid, 'compile(%u) \n', nChains);

initfileN = size(initStructs,2);
for i=1:initfileN
  initFileName = fullfileKPM(workingDir, ['init_', num2str(i) '.txt']);
  dataGen(initStructs(i), initFileName)
  fprintf(fid, 'inits (%u, ''%s'')\n', i, initFileName);
end

if 0 % blocking
  fprintf(fid, 'blockfe(1)\n'); % block fixed effects (for GLMs)
end
fprintf(fid, 'gen.inits() \n');
fprintf(fid, 'update(%u)\n', nBurnin);

%setting params to monitor
if isempty(monitorParams)
  fprintf(fid, 'set (*)\n');
else
  for i=1:length(monitorParams)
    fprintf(fid, 'set (%s)\n', strrep(monitorParams{i}, '_', '.'));
  end
  %fprintf(fid, 'set (%s)\n','deviance');
end


%fprintf(fid, 'dic.set()\n');
fprintf(fid, 'thin.updater(%u)\n', thin);
fprintf(fid, 'update(%u)\n', nSamples);
fprintf(fid, 'coda(*, ''%s'')\n',  codaFile);
fprintf(fid, 'stats(*)\n');


%fprintf(fid, 'dic.stats()\n');
fprintf(fid, 'history(*, ''%s'')\n', his_filename);
fprintf (fid, 'save (''%s'')\n',  log_filename);

if (view == 0)
  fprintf(fid, 'quit() \n');
end
fclose(fid);

%calling WinBUGS and passing the script to it
str = ['"', Bugdir, '\Winbugs14.exe" /PAR script.txt'];
dos(str); % DOS wants a\b

% passing the coda files to bugs2mat to convert files to structs
codaIndex = [codaFile, 'Index.txt'];
for i=1:nChains
  codaF = [codaFile, num2str(i), '.txt'];
  S = bugs2mat(codaIndex, codaF);
  structArray(i) = S;
end

samples = structsToArrays(structArray);
stats = computeStats(samples);

if nChains == 1
   disp('EPSR not calculated (only one chain)');
end

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

function dataGen(dataStruct, fileName)
% This is a helper function to generate data or init files for winBUGS
% Inputs:
%   fileName: name of the text file containing initial values. for each
%             chain we'll fileName_i where 'i' is the chain number,
%   dataStruct: is a Struct with name of params(consistant in the same
%               order with paramList) are fields and intial values are functions

if nargin<2
   error(['This function needs two arguments']);
end

fieldNames = fieldnames(dataStruct);
Nparam = size(fieldNames, 1);

%fileName = [fileName, '.txt'];
fid = fopen(fileName, 'w');
if fid == -1
  error(['Cannot open ', fileName ]);
end

fprintf(fid,'list(');
for i=1:Nparam
  fn = fieldNames(i);
  fval = fn{1};
  val = getfield(dataStruct, fval);
  [sfield1, sfield2]= size(val);

  msfield = max(sfield1, sfield2);
  newfval = strrep(fval, '_', '.');
  if ((sfield1 == 1) && (sfield2 == 1))  % if the field is a singleton
    fprintf(fid, '%s=%G',newfval, val);

  %
  % One-D array:
  %   beta = c(6, 6, ...)
  % 
  % 2-D or more:
  %   Y=structure(
  %     .Data = c(1, 2, ...), .Dim = c(30,5))
  %
  elseif ((length(size(val)) == 2) && ((sfield1 == 1) || (sfield2 == 1)))
    fprintf(fid, '%s=c(',newfval);
    for j=1:msfield
      if (isnan(val(j)))
        fprintf(fid,'NA');
      else
        % format for winbugs
        fprintf(fid,wb_strval(val(j)));
      end
      if (j<msfield)
        fprintf(fid, ', ');
      else
        fprintf(fid, ')');
      end
    end
  else
    % non-trivial 2-D or more array
    valsize    = size(val);
    alldatalen = prod(valsize);
    alldata = reshape(val', [1, alldatalen]);
    fprintf(fid, '%s=structure(.Data=c(', newfval);
    for j=1:alldatalen
      if (isnan(alldata(j)))
        fprintf(fid,'NA');
      else
        % format for winbugs
        fprintf(fid,wb_strval(alldata(j)));
      end
      if (j < alldatalen)
        fprintf(fid,',');
      else
        fprintf(fid,'), .Dim=c(', alldata(j));
      end
    end

    for j=1:length(valsize)
      if (j < length(valsize))
        fprintf(fid, '%G,', valsize(j));
      else
        fprintf(fid, '%G))', valsize(j));
      end
    end
  end
  if (i<Nparam)
    fprintf(fid, ', ');
  else
    fprintf(fid, ')\n');
  end
end
fclose(fid);

%%%%%%%%

function s = wb_strval(v)
% Converts numeric value to a string that is acceptable by winbugs.
% This is most problematic for exponent values which must have at least 1
% decimal and two digits for the exponent. So 1.0E+01 rather than 1E+001
% Note that only Matlab on PC does 3 digits for exponent.
s = sprintf('%G', v);
if strfind(s, 'E')
   if length(strfind(s, '.')) == 0
      s = strrep(s, 'E', '.0E');
   end
   s = strrep(s, 'E+0', 'E+');
   s = strrep(s, 'E-0', 'E-');
end

%%%%%%%%

function f = fullfileKPM(varargin)
% fullfileKPM Concatenate strings with file separator, then convert it to a/b/c
% function f = fullfileKPM(varargin)

f = fullfile(varargin{:});
f = strrep(f, '\', '/');

%%%%%%%%

function A = structsToArrays(S)
% Suppose S is this struct array
%
% S(c).X1(s)
% S(c).X2(s,i)
% S(c).X3(s,i,j)
%
% where s=1:N in all cases 
%
% Then we return
% A.X1(c,s)
% A.X2(c,s,i)
% A.X3(c,s,i,j)

C = length(S);
fld = fieldnames(S);
A = [];
for fi=1:length(fld)