Skip to content

Commit 6068f42

Browse files
NilsWinterNilsWinter
committed
added one group estimation
1 parent f0e9a76 commit 6068f42

File tree

7 files changed

+273
-2
lines changed

7 files changed

+273
-2
lines changed

‎README.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -162,6 +162,15 @@ These are examples of the figures that can be created to show the posterior dist
162162
* mbe_2gr_summary.m
163163
> - Computes summary statistics for all parameters of a 2 group comparison.This will only work for a mcmc chain with parameters mu1,mu2,sigma1,sigma2 and nu.
164164
165+
* mbe_1gr_example.m
166+
> - This is an example script for a one group Bayes estimation.
167+
168+
* mbe_1gr_plots.m
169+
> - Makes histogram of data with superimposed posterior prediction check and plots posterior distribution of monitored parameters.
170+
171+
* mbe_1gr_summary.m
172+
> - Computes summary statistics for all parameters.This will only work for a mcmc chain with parameters mu1,sigma1 and nu.
173+
165174
#### MCMC Diagnostics
166175
* mbe_diagMCMC.m
167176
> - Plots autocorrelation, parameter trace, shrink factor and parameter density.

‎exampleFigs/plots_for_1gr_example.tif‎

2.98 MB
Binary file not shown.

‎mbe_1gr_example.m‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
%% Matlab Toolbox for Bayesian Estimation - One Group Example
2+
% This is an example script for a one group estimation.
3+
4+
% Largely based on R code introduced in the following paper:
5+
% Kruschke, J.K., Bayesian Estimation supersedes the t-test.
6+
% Journal of Experimental Psychology: General, Vol 142(2), May 2013, 573-603.
7+
% see http://www.indiana.edu/~kruschke/BEST/ for R code
8+
% Nils Winter ([email protected])
9+
% Johann-Wolfgang-Goethe University, Frankfurt
10+
% Created: 2016-04-25
11+
% Version: v1.0 (2016-04-25)
12+
%-------------------------------------------------------------------------
13+
clear;clc;
14+
15+
%% Load some data
16+
% EXAMPLE DATA (see Kruschke, 2013)
17+
% Run this script one time with the small sample.
18+
y = [101,100,102,104,102,97,105,105,98,101,100,123,105,103,100,95,102,106,...
19+
109,102,82,102,100,102,102,101,102,102,103,103,97,97,103,101,97,104,...
20+
96,103,124,101,101,100,101,101,104,100,101];
21+
nTotal = length(y);
22+
23+
% % Then run it again but use the n=500 sample below to see the change
24+
% % of the parameters and HDI.
25+
% y = trnd(2,500,1);
26+
% y = (y-mean(y))/sqrt(mean((y-mean(y)).^2))*6+102;
27+
% nTotal = length(y);
28+
29+
%% Specify prior constants, shape and rate for gamma distribution
30+
% See Kruschke (2013) for further description
31+
muM = mean(y);
32+
muP = 0.000001 * 1/std(y)^2;
33+
sigmaLow = std(y)/1000;
34+
sigmaHigh = std(y)*1000;
35+
36+
% Save prior constants in a structure for later use with matjags
37+
dataList = struct('y',y,'nTotal',nTotal,...
38+
'muM',muM,'muP',muP,'sigmaLow',sigmaLow,'sigmaHigh',sigmaHigh);
39+
40+
%% Specify MCMC properties
41+
% Number of MCMC steps that are saved for EACH chain
42+
% This is different to Rjags, where you would define the number of
43+
% steps to be saved for all chains together (in this example 12000)
44+
numSavedSteps = 30000;
45+
46+
% Number of separate MCMC chains
47+
nChains = 3;
48+
49+
% Number of steps that are thinned, matjags will only keep every nth
50+
% step. This does not affect the number of saved steps. I.e. in order
51+
% to compute 10000 saved steps, matjags/JAGS will compute 50000 steps
52+
% If memory isn't an issue, Kruschke recommends to use longer chains
53+
% and no thinning at all.
54+
thinSteps = 1;
55+
56+
% Number of burn-in samples
57+
burnInSteps = 1000;
58+
59+
% The parameters that are to be monitored
60+
parameters = {'mu','sigma','nu'};
61+
62+
%% Initialize the chain
63+
% Initial values of MCMC chains based on data:
64+
mu = mean(y);
65+
sigma = std(y);
66+
% Regarding initial values: (1) sigma will tend to be too big if
67+
% the data have outliers, and (2) nu starts at 5 as a moderate value. These
68+
% initial values keep the burn-in period moderate.
69+
70+
% Set initial values for latent variable in each chain
71+
for i=1:nChains
72+
initsList(i) = struct('mu', mu, 'sigma',sigma,'nu',5);
73+
end
74+
75+
%% Specify the JAGS model
76+
% This will write a JAGS model to a text file
77+
% You can also write the JAGS model directly to a text file
78+
79+
modelString = [' model {\n',...
80+
' for ( i in 1:nTotal ) {\n',...
81+
' y[i] ~ dt( mu , tau, nu )\n',...
82+
' }\n',...
83+
' mu ~ dnorm( muM , muP ) \n',...
84+
' tau <- 1/pow(sigma , 2)\n',...
85+
' sigma ~ dunif( sigmaLow , sigmaHigh )\n',...
86+
' nu ~ dexp( 1/30 )\n'...
87+
'}'];
88+
fileID = fopen('mbe_1gr_example.txt','wt');
89+
fprintf(fileID,modelString);
90+
fclose(fileID);
91+
model = fullfile(pwd,'mbe_1gr_example.txt');
92+
93+
%% Run the chains using matjags and JAGS
94+
% In case you have the Parallel Computing Toolbox, use ('doParallel',1)
95+
[~, ~, mcmcChain] = matjags(...
96+
dataList,...
97+
model,...
98+
initsList,...
99+
'monitorparams', parameters,...
100+
'nChains', nChains,...
101+
'nBurnin', burnInSteps,...
102+
'thin', thinSteps,...
103+
'verbosity',1,...
104+
'nSamples',numSavedSteps);
105+
106+
%% Restructure the output
107+
% This transforms the output of matjags into the format that mbe is
108+
% using
109+
mcmcChain = mbe_restructChains(mcmcChain);
110+
111+
%% Examine the chains
112+
mbe_diagMCMC(mcmcChain);
113+
114+
%% Examine the results
115+
% At this point, we want to use all the chains at once, so we
116+
% need to concatenate the individual chains to one long chain first
117+
mcmcChain = mbe_concChains(mcmcChain);
118+
% Get summary and posterior plots; comparisonValue = 100
119+
summary = mbe_1gr_summary(mcmcChain,100);
120+
% Data has to be in a cell array and the vectors have to be column vectors
121+
data = {y'};
122+
mbe_1gr_plots(data,mcmcChain,100);

‎mbe_1gr_example.txt‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
model {
2+
for ( i in 1:nTotal ) {
3+
y[i] ~ dt( mu , tau, nu )
4+
}
5+
mu ~ dnorm( muM , muP )
6+
tau <- 1/pow(sigma , 2)
7+
sigma ~ dunif( sigmaLow , sigmaHigh )
8+
nu ~ dexp( 1/30 )
9+
}

‎mbe_1gr_plots.m‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
function mbe_1gr_plots(y, mcmcChain, compVal, varargin)
2+
% mbe_makePlots
3+
% Make histogram of data with superimposed posterior prediction check
4+
% and plots posterior distribution of monitored parameters.
5+
%
6+
% INPUT:
7+
% y
8+
% cell array containing vectors for y1
9+
% mcmcChain
10+
% structure with one MCMC-chain, should contain all monitored parameters
11+
% compVal
12+
% comparison value for computing of effect size
13+
%
14+
% Specify the following name/value pairs for additional plot options:
15+
% Parameter Value
16+
% 'plotPairs' show correlation plot of parameters ([1],0)
17+
%
18+
%
19+
% EXAMPLE:
20+
21+
% Largely based on R code introduced in the following paper:
22+
% Kruschke, J.K., Bayesian Estimation supersedes the t-test.
23+
% Journal of Experimental Psychology: General, Vol 142(2), May 2013, 573-603.
24+
% see http://www.indiana.edu/~kruschke/BEST/ for R code
25+
% Nils Winter ([email protected])
26+
% Johann-Wolfgang-Goethe University, Frankfurt
27+
% Created: 2016-04-25
28+
% Version: v1.00 (2016-04-25)
29+
%-------------------------------------------------------------------------
30+
31+
% -----------------------------------------------------------------
32+
% Get input
33+
% -----------------------------------------------------------------
34+
p = inputParser;
35+
defaultPlotPairs = 1;
36+
addOptional(p,'plotPairs',defaultPlotPairs);
37+
parse(p,varargin{:});
38+
plotPairs = p.Results.plotPairs;
39+
40+
% Get parameter names
41+
names = fieldnames(mcmcChain);
42+
43+
%% -----------------------------------------------------------------
44+
% Plot correlations between parameters
45+
%-----------------------------------------------------------------
46+
if plotPairs
47+
mbe_plotPairs(mcmcChain,1000)
48+
end
49+
50+
%% -----------------------------------------------------------------
51+
% Plot data y and smattering of posterior predictive curves:
52+
%-----------------------------------------------------------------
53+
nu = mcmcChain.(names{3});
54+
mu = mcmcChain.(names{1});
55+
sigma = mcmcChain.(names{2});
56+
figure('Color','w','NumberTitle','Off','Position',[100,50,800,600]);
57+
subplot(3,2,[2 4]);
58+
mbe_plotData(y,nu,mu,sigma);
59+
60+
%% -----------------------------------------------------------------
61+
% Plot posterior distribution of parameter nu:
62+
%-----------------------------------------------------------------
63+
subplot(3,2,6);
64+
mbe_plotPost(log10(nu),'credMass',0.95,'xlab','log10(\nu)','PlotTitle','Normality');
65+
66+
%-----------------------------------------------------------------
67+
% Plot posterior distribution of parameters mu:
68+
%-----------------------------------------------------------------
69+
xLim(1) = min(mu);
70+
xLim(2) = max(mu);
71+
72+
subplot(3,2,1);
73+
mbe_plotPost(mu,'xlab','\mu','xlim',xLim,'Plottitle','Mean');
74+
75+
%-----------------------------------------------------------------
76+
% Plot posterior distribution of param's sigma1, sigma2, and their difference:
77+
%-----------------------------------------------------------------
78+
xLim(1) = min(sigma);
79+
xLim(2) = max(sigma);
80+
subplot(3,2,3);
81+
mbe_plotPost(sigma,'xlab','\sigma','xlim',xLim,'PlotTitle','Std. Dev.');
82+
83+
%-----------------------------------------------------------------
84+
% Plot of estimated effect size. Effect size is d-sub-a from
85+
%-----------------------------------------------------------------
86+
% Macmillan & Creelman, 1991; Simpson & Fitter, 1973; Swets, 1986a, 1986b.
87+
effectSize = (mu - compVal) ./ sigma;
88+
subplot(3,2,5);
89+
str = ['(\mu-' num2str(compVal) ')/\sigma'];
90+
mbe_plotPost(effectSize,'rope',[-0.1,0.1],'xlab',str,'PlotTitle','Effect Size');
91+
92+
93+
end

‎mbe_1gr_summary.m‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function summary = mbe_1gr_summary(mcmcChain,compVal)
2+
%% mbe_1gr_summary
3+
% Computes summary statistics for all parameters of a one group estimation.
4+
% This will only work for a mcmc chain with parameters mu1,sigma1,
5+
% and nu.
6+
%
7+
% INPUT:
8+
% mcmcChain
9+
% structure with fields for mu, sigma, nu
10+
% compVal
11+
% comparison value. Needed to compute effect size.
12+
%
13+
% OUTPUT:
14+
% summary
15+
% outputs structure containing mu1, sigma1,
16+
% nu, nuLog10 and effectSize
17+
%
18+
% EXAMPLE:
19+
% summary = mbe_1gr_summary(mcmcChain);
20+
21+
% Largely based on R code introduced in the following paper:
22+
% Kruschke, J.K., Bayesian Estimation supersedes the t-test.
23+
% Journal of Experimental Psychology: General, Vol 142(2), May 2013, 573-603.
24+
% see http://www.indiana.edu/~kruschke/BEST/ for R code
25+
% Nils Winter ([email protected])
26+
% Johann-Wolfgang-Goethe University, Frankfurt
27+
% Created: 2016-04-25
28+
% Version: v1.0 (2016-04-25)
29+
%-------------------------------------------------------------------------
30+
summary.mu = mbe_summary(mcmcChain.mu1);
31+
summary.sigma = mbe_summary(mcmcChain.sigma1);
32+
summary.nu = mbe_summary(mcmcChain.nu1);
33+
summary.nuLog10 = mbe_summary(log10(mcmcChain.nu1));
34+
effSzChain = (mcmcChain.mu1 - compVal)./mcmcChain.sigma1;
35+
summary.effSz = mbe_summary(effSzChain,0);
36+
end
37+

‎mbe_plotData.m‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,9 @@ function mbe_plotData(y,nu,mu1,sigma1,mu2,sigma2)
4747
end
4848

4949
%% Plot data for y1
50-
figure('Color','w','NumberTitle','Off','Position',[100,50,800,600]);
51-
subplot(size(y,2),1,1);
50+
if exist('y2','var')
51+
subplot(2,1,1);
52+
end
5253
% Select thinned steps in chain for plotting of posterior predictive curves:
5354
chainLength = size(mu1,1);
5455
nCurvesToPlot = 30;

0 commit comments

Comments
 (0)