Dear All:

---
UPDATE: I have made some adaptations to the original matbugs.m script and it now works on Wine with the Rate_1.m example from the wonderful book below. I have emailed the authors and am hoping that my small contributions (marked with my name in the file) might help others. Thank you so much.

My original message included below.

Misha

---

Thank you so much for all your wonderful input on WinBUGS14. I am currently using this great product with Wine (www.winehq.org) on Ubuntu Hardy Heron 8.04. I will try to summarize all your prior responses once I have been able to implement them successfully for our application.

I have found the MatBUGS package found here:
http://people.cs.ubc.ca/~murphyk/Software/MATBUGS/matbugs.html

to be quite helpful, as relayed by the free book linked from this page:
http://web.mit.edu/yarden/www/bayes.html

I was wondering if there was a version of MATBUGS available for WinBUGS14 run via Wine?

Thank you so much

Sincerely yours! Happy Thanksgiving!
Misha Koshelev
-- 
Misha Koshelev
MD/PhD Student

Human Neuroimaging Laboratory
One Baylor Plaza
S104
Baylor College of Medicine
Houston, TX 77030


-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list



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]
% 'openBugs' - set to 1 to use openBugs file format [0]
% 'Bugdir'   - location of winbugs executable
%               Default is 'C:/Program Files/WinBUGS14' if not openBugs
%               Default is 'C:/Program Files/OpenBUGS' if OpenBugs.
% 'workingDir' - directory to store temporary data/init/coda files [pwd/tmp]
%
% 'DICstatus' - takes value 1 to set the DIC tool and 0 otherwise
% 'refreshrate' - sets the refresh rate for the updater. Default is 100. Values
%               of 10 prevent the computer from hanging too much if the model
%               is very slow to run.
%
% 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, stats.DIC (Andrew Jackson added).
% Each field may have a different shape:
%    stats.mean.theta
%    stats.mean.theta(i)
%    stats.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.
%
% DIC reads the actual DIC values generated by WinBUGS. This is different
%   to Andrew Gelman's R2WinBUGS which calculates these stats by monitoring
%   the deviance. Problem with his method is that he must estimate pD by
%   var(deviance)/2. This can be quite different to the winbugs value.
%   matbugs reads the DIC values from the winbugs log file directly. Fields
%   are generated for each of the variables listed including the total value.
%   e.g. stats.DIC.total returns the fields Dbar, Dhat, pD, DIC as per
%   winbugs. So stats.DIC.total.DIC returns the actual DIC value you
%   probably want.
%
% Example
%
%     [S,stats] = matbugs(data, 'schools_model.txt', 'nSamples', 10000, ...
%                          'monitorParams', {'mu_theta', 'theta'}, ...
%                          'init', initStruct, 'nchains', 3, 'DICstatus',1)
%
% Written by Maryam Mahdaviani ([log in to unmask])
% and Kevin Murphy ([log in to unmask]), August 2005
% Modified for OpenBUGS by Tomo Eguchi 23 March 2006
% Modified for DIC by Andrew Jackson ([log in to unmask]), March 2006 
% Modified for winbugs filename by  Sohrab Shah 28 Nov 2006
% Bug fix 7 Mar 08 Woojae Kim clear valTransp

[openBUGS, junk] =  process_options(varargin, 'openBUGS', 0);

if openBUGS
  Bugdir = 'C:/Program Files/OpenBUGS';
else
  Bugdir = 'C:/Program Files/WinBUGS14';
end

% Misha Koshelev Nov 26 2009
% if on unix, assume wine
% if wine, z: = /; otherwise leave blank
if isunix
    wine=1;
    pathPrefix='z:';
else
    wine=0;
    pathPrefix='';
end

[initStructs, Bugdir, nChains, view, workingDir, nBurnin, nSamples, ...
 monitorParams, thin, blocking, refreshrate, DICstatus, openBUGS, junk] =  ...
   process_options(...
	varargin, ...
	'init', {}, ...
	'Bugdir', Bugdir, ...
	'nChains', 3, ...
	'view', 0, ...
	'workingDir', fullfile(pwd,'tmp'), ...
	'nBurnin', 1000, ...
	'nSamples', 5000, ...
	'monitorParams', {}, ...
	'thin', 1, ...
	'blocking', 1, ...
       'refreshrate',100,...
       'DICstatus',0, ...
       'openBUGS', 0);

     % allow user to not specify initial values
if 0 % 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');

% fixed for wine Misha Koshelev Nov 26 2009
if wine
    scriptFile = [Bugdir,filesep,'script.txt'];
else
    scriptFile = [Bugdir,'\','script.txt'];
end

%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)
if openBUGS
  fprintf(fid, 'modelCheck(''%s'')\n',[pathPrefix,bugsModel]); % Misha Koshelev Nov 26th 2009
else
  fprintf(fid, 'check(''%s'')\n',[pathPrefix,bugsModel]); % Misha Koshelev Nov 26th 2009
end


if ~isempty(dataStruct)
  dataFile = fullfileKPM(workingDir, 'data.txt');
  dataGen(dataStruct, dataFile);
  if openBUGS
    fprintf(fid, 'modelData(''%s'')\n', [pathPrefix,dataFile]); % Misha Koshelev Nov 26th 2009
  else
    fprintf(fid, 'data(''%s'')\n', [pathPrefix,dataFile]); % Misha Koshelev Nov 26th 2009
  end
end
%fprintf(fid, '.txt'')\n');

if openBUGS
  fprintf(fid, 'modelCompile(%u) \n', nChains);
else
  fprintf(fid, 'compile(%u) \n', nChains);
end    


initfileN = size(initStructs,2);
for i=1:initfileN
  initFileName = fullfileKPM(workingDir, ['init_', num2str(i) '.txt']);
  dataGen(initStructs(i), initFileName)
  if openBUGS,
      fprintf(fid, 'modelInits(''%s'', %u)\n', [pathPrefix,initFileName], i); % Misha Koshelev Nov 26th 2009
 else
     fprintf(fid, 'inits (%u, ''%s'')\n', i, [pathPrefix,initFileName]); % Misha Koshelev Nov 26th 2009
 end
end


if 0 
  % block fixed effects (for GLMs)
  %    Don't know how to make this work in winbugs
  % Not available in openbugs
  fprintf(fid, 'blockfe(1)\n'); 
end

fprintf(fid, 'refresh(%u) \n', refreshrate); % Andrew Jackson - line added

if openBUGS
  fprintf(fid, 'modelGenInits() \n');
  fprintf(fid, 'modelUpdate(%u, TRUE)\n', nBurnin);
else
  fprintf(fid, 'gen.inits() \n');
  fprintf(fid, 'update(%u)\n', nBurnin);
end

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

if DICstatus; fprintf(fid, 'dic.set()\n'); end % Andrew Jackson - line added

if openBUGS
  fprintf(fid, 'samplesThin(%u)\n', thin);
  fprintf(fid, 'modelUpdate(%u)\n', nSamples);
  fprintf(fid, 'samplesCoda("*", ''%s'')\n',  [pathPrefix,codaFile]); % Misha Koshelev Nov 26th 2009
  fprintf(fid, 'samplesStats("*")\n');
  fprintf(fid, 'samplesDensity("*")\n');
  fprintf(fid, 'samplesHistory("*")\n');
else
  fprintf(fid, 'thin.updater(%u)\n', thin);
  fprintf(fid, 'update(%u)\n', nSamples);
  fprintf(fid, 'coda(*, ''%s'')\n',  codaFile);
  fprintf(fid, 'stats(*)\n');
end

% Andrew Jackson - line added, includes marker text to denote the end of
% the DIC stats see below
if DICstatus; fprintf(fid, 'dic.stats()\n #endDIC'); end

if openBUGS
  fprintf(fid, 'samplesHistory("*", ''%s'')\n', his_filename);
  fprintf (fid, 'modelSaveLog(''%s'')\n',  log_filename);
else
  fprintf(fid, 'history(*, ''%s'')\n', his_filename);
  fprintf (fid, 'save (''%s'')\n',  log_filename);
end

if (view == 0)
  if openBUGS
    %fprintf(fid, 'modelQuit() \n');
    % Bug fix by Brani Vidakovic
    fprintf(fid, 'modelQuit("y")\n');
  else
    fprintf(fid, 'quit() \n');
  end
end
fclose(fid);


%calling WinBUGS and passing the script to it
% File name fix by Sohrab Shah 28 Nov 2006
if openBUGS
 f = fullfile(Bugdir, 'winbugs.exe');
 %str = ['"', Bugdir, '\winbugs.exe" /PAR script.txt'];
else
 % fix for wine Misha Koshelev 26 Nov 2009
 if wine
     f = fullfile(Bugdir, 'WinBUGS14.exe');
 else
     f = fullfile(Bugdir, 'Winbugs14.exe');
 end
 %str = ['"', Bugdir, '\Winbugs14.exe" /PAR script.txt'];
end
str = ['"',f,'" /PAR script.txt'];
dos(str); % DOS wants a\b


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

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

if DICstatus;
  DICstats = getDICstats(workingDir);
  stats.DIC = DICstats; 
end

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]);
    %alldata = alldata(:)';
    
    %Truccolo-Filho, Wilson <[log in to unmask]>
   if length(valsize)<3
     alldata = reshape(val', [1, alldatalen]);
   elseif length(valsize)==3
     clear valTransp
     for j=1:valsize(3)
       valTransp(j,:,:)=val(:,:,j)';%need a new variable, since val might be rectangular
     end
     alldata=valTransp(:)';
   else
     ['Error: 4D and higher dimensional arrays not accepted']
     return
   end
    
    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)
  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)
%
% 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  [log in to unmask]
  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

%%%%%%%%%%%%
% Andrew Jackson - function getDICstats added [log in to unmask]
% Used to retrieve the DIC statistics from the log file
% if matbugs input DICstatus = 1
% This code is probably a bit clunky but it does the job.
% Note that the values read from the log file have 3 decimal places but
% matlab decides to give them 4 - with a zero tagged on the end. The
% precision is obviously only to 3 decimal places. sscanf.m does not seem
% to recognise the field-width and precision codes that sprintf.m does.
function DICstats = getDICstats(workingDir)
DICstats = [];
FIDlog = fopen([workingDir '\log.txt'],'r');
ct = 0;
test = 0;
endloop = 0;
while 1
    
    tline = fgets(FIDlog);
    
    if tline == -1; break; end
    if endloop; break; end
    
    if strfind(tline,'dic.set cannot be executed'); 
        DICstats.error = 'DIC monitor could not be set by WinBUGS';
    end
    
    if size(tline,2)>6
        % The string 'total' in the log file denotes the end of the DIC
        % stats so the loop can be ended in the next iteration.
        if strcmp(tline(1:5),'total'); endloop = 1; end;
    end

    if size(tline,2)>2
        % locate the DIC string identifier in the log file
        if strcmp(tline(1:3),'DIC'); test = 1; end
    end
    
    if test 
        ct=ct+1;
        % DIC results are located 3 lines after the DIC string identifier
        % in the log file.
        if ct >= 4
            A = sscanf(tline,'%*s %f %f %f %f');
            S = sscanf(tline, '%s %*f %*f %*f %*f');
            % Cheng-Ta, Yang suggested this change
            %DICstats = setfield(DICstats,S,'Dbar',A(1));
            %%DICstats = setfield(DICstats,S,'Dhat',A(2));
            %DICstats = setfield(DICstats,S,'pD',A(3));
            %DICstats = setfield(DICstats,S,'DIC',A(4));
            DICstats.S.Dbar = A(1);
            DICstats.S.Dhat = A(2);
            DICstats.S.pD = A(3);
            DICstats.S.DIC = A(4);
            DICstats.S.Dhat = A(2);
        end
    end
end

fclose(FIDlog)


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

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) [log in to unmask], 2000
% 2003-01-14 [log in to unmask] - 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
% GNU GPL

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


-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list



% Data
k=5;n=10;
% WinBUGS Parameters
nchains=2; % How Many Chains?
nburnin=0; % How Many Burn-in Samples?
nsamples=2e4; %How Many Recorded Samples?

% Assign Matlab Variables to the Observed WinBUGS Nodes
datastruct = struct('k',k,'n',n);
% Initialize Unobserved Variables
start.theta= [0.1 0.9];
for i=1:nchains
S.theta = start.theta(i); % An Intial Value for the Success Rate
init0(i) = S;
end
% Use WinBUGS to Sample
[samples, stats] = matbugs(datastruct, ...
fullfile(pwd, 'Rate_1.txt'), ...
'init', init0, ...
'nChains', nchains, ...
'view', 1, 'nburnin', nburnin, 'nsamples', nsamples, ...
'thin', 1, 'DICstatus', 0, 'refreshrate',100, ...
'monitorParams', {'theta'}, ...
'Bugdir', '/home/misha/.wine/drive_c/Program Files/WinBUGS14');