Beta_MPT_Model

Both the Model and MatLab Code displayed here can be found in the attachements at the bottom of this page.

To run this 
  • You would need to first download the free software WinBUGS at: 
    http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml

  • Secondly, I am using MatLab so to give MatLab the ability to interact with WinBUGS, download the freely available matbugs.m function and put it in your Matlab working directory. You can download matbugs.m directly from http://www.cs.ubc.ca/~murphyk/Software/MATBUGS/ matbugs.html


# Save this model in a text file as Beta_MPT_Pair_Clustering.txt and you will be able to run the following code found below the model in #MatLab/WinBUGS

model {    
    # Following loop defines the base MPT model
 
for (n in 1:subjs) {   
    # Loop iterates through, subjs, subjects
    # theta[s,n] defines the base parameters where s indexes the parameters, and n the subjects
    # p[n,k] defines the category probabilities where, n indexes the subject and k the category
 
    p[n,1] <- theta[1,n]*theta[2,n]
    p[n,2] <- (1-theta[1,n])*pow(theta[3,n],2)
    p[n,3] <- (1-theta[1,n])*(2*theta[3,n])*(1-theta[3,n])
    p[n,4] <- theta[1,n]*(1-theta[2,n])+(1-theta[1,n])*pow((1-theta[3,n]),2)
 
    # Designates the data follows a multinomial distribution given the category probabilities response[n,k] is the data matrix
    response[n,1:4] ~ dmulti(p[n,1:4],items)
}
 
# Following loop defines the hierarchical distribution 
for (s in 1:3) {
    for (n in 1:subjs) {
    theta[s,n] ~ dbeta(alph[s],bet[s])I(.001,.999)   # Designates that the base MPT parameters are beta distributed
    }
}
 
# Priors for each parameter
for (s in 1:3) {
 # Sets up the priors to be close to flat on the mean and standard deviation of the beta distribution
    zero[s] <- 0
    zero[s] ~ dpois(phi[s])
    phi[s] <- -log(1/pow(alph[s]+bet[s],5/2)) 
    
    alph[s] ~ dunif(.01,5000)
    bet[s] ~ dunif(.01,5000)
    
    mnb[s] <- alph[s]/(alph[s]+bet[s]) # Used to monitor the mean
    varp[s] <- sqrt(alph[s]*bet[s]/(pow(alph[s]+bet[s],2)*(alph[s]+bet[s]+1)))
    }
}

_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+_+



% This is what you would add in Matlab after you had installed MatBUGS to run WinBUGS. The data set that you see here is Trial 1 #from Organic Alcoholics. 

response=[14   3   1   2
12  3   1   4
18  0   1   1
15  3   0   2
7   1   10  2
3   6   11  0
8   4   3   5
17  1   1   1
13  4   3   0
11  6   1   2
16  1   2   1
10  1   3   6
7   13  0   0
8   4   3   5
16  1   1   2
5   4   7   4
15  0   5   0
6   3   6   5
17  2   0   1
17  1   0   2
8   3   6   3];
 
[subjs notused] =size(response);
 
items=20; % The number of items used in the experiment. A row would add up to 20

 
nchains = 1; % number of chains
nburnin = 20000; % number of burn-in samples
nsamples = 25000; % number of samples
 
% Data to Supply to WinBugs
datastruct = struct('subjs',subjs,'response',response,'items',items);
 
% Initial Values to Supply to WinBugs
for i=1:nchains
 S.alph = 1/2;
 S.bet = 1/2;
 S.theta = 1/2;
 init0(i) = S;
end
 
 [samples, stats] = matbugs(datastruct, ...
 fullfile(pwd, 'Beta_MPT_Pair_Clustering.txt'), ...
 'init', init0, ...
 'nChains', nchains, ...
 'view', 1, 'nburnin', nburnin, 'nsamples', nsamples, ...
 'thin', 1, 'DICstatus', 0, 'refreshrate',100, ...
 'monitorParams', {'alph','bet','varp','mnb','lambda'}, ...
 'Bugdir', 'C:/Program Files/WinBUGS14');

ċ
Beta_MPT_Pair_Clustering.txt
(2k)
Gregory Alexander,
Jul 19, 2011, 6:07 AM
ċ
Beta_MPT_Pair_Clustering_Model.m
(1k)
Gregory Alexander,
Jul 19, 2011, 6:14 AM
Comments