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 (maryam@cs.ubc.ca)% and Kevin Murphy (murphyk@cs.ubc.ca), 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)]);endif ~exist(workingDir, 'dir') mkdir(pwd, 'tmp');endlog_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\bcodaFile = 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);endif 0 % blocking fprintf(fid, 'blockfe(1)\n'); % block fixed effects (for GLMs)endfprintf(fid, 'gen.inits() \n');fprintf(fid, 'update(%u)\n', nBurnin);%setting params to monitorif 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');endfclose(fid);%calling WinBUGS and passing the script to itstr = ['"', Bugdir, '\Winbugs14.exe" /PAR script.txt'];dos(str); % DOS wants a\b% passing the coda files to bugs2mat to convert files to structscodaIndex = [codaFile, 'Index.txt'];for i=1:nChains codaF = [codaFile, num2str(i), '.txt']; S = bugs2mat(codaIndex, codaF); structArray(i) = S;endsamples = 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 functionsif nargin error(['This function needs two arguments']);endfieldNames = fieldnames(dataStruct);Nparam = size(fieldNames, 1);%fileName = [fileName, '.txt'];fid = fopen(fileName, 'w');if fid == -1 error(['Cannot open ', fileName ]);endfprintf(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 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 fprintf(fid, ', '); else fprintf(fid, ')\n'); endendfclose(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)
443
评论
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
查看更多评论
添加红包
