Navigation

Personal Bayesian Simple Reaction Time Model

WebPPL-Code

console.log("==============================================================================")
console.log("PCM20180306_PersonalBayesianSimpleReactionTimeModel *** version 2018/03/27 ***")
console.log("see also Simple Reaction Time, Example 9, Card, Moran & Newell, 1983, p.66    ")
console.log("see also www.humanbenchmark.com/tests/reactiontime/statistics ")
console.log("==============================================================================")
/**
 * @author - Claus Moebus   <claus.moebus@uol.de>
 */
//------------------------------------------------------------------------------------------
/**
 * @variable {number} startTime - used in method 'runtime' to compute runtime in sec and min
 */
var startTime = Date.now()
//------------------------------------------------------------------------------------------
console.log("Input parameter:")
/**
 * @variable {integer} nTrials - no of efficient samples (incl. burnout) in MCMC-sampling
 */
var nTrials = 20E3
console.log("nTrials = ", nTrials)
//------------------------------------------------------------------------------------------
/**
 * @variable {integer} nSigma - no of standard deviations between mean and 'slow', 'fast'
 *                              interval boundaries
 */
var nSigma = 1
console.log("nSigma = ", nSigma)
//------------------------------------------------------------------------------------------
/**
 * @variable {integer} myBurnPeriod - length of burnin period in MCMC process
 */
var myBurnPeriod = nTrials * 0.10
console.log("length of burn-in period = ", myBurnPeriod)
//------------------------------------------------------------------------------------------
/**
 * @variable {integer} myLag - only every myLag-th sample will be retained during MCMC
 */
var myLag = 10
console.log("length of lag = ", myLag)
//------------------------------------------------------------------------------------------
/**
 * @variable {array} data - author's reaction times in an experiment found here
 *                          www.humanbenchmark.com/tests/reactiontime/
 *                          visited March 2018
 */
var data =                                        
    [458, 292, 228, 403, 271, 420, 350, 235, 260, 306]
console.log("response time data = [", data, "]")
//------------------------------------------------------------------------------------------
/**
 * @object hyperParmTauT_a - contains parameters for the uniform prior of tauT_a
 * @property {number} a - value is the lower boundary of the uniform density
 * @property {number} b - value is the upper boundary of the uniform density
 */
var hyperParmTauT_a = {a:5.00, b: 50.00}
console.log("hyperParmTauT_a.a = ", hyperParmTauT_a.a,
            ", hyperParmTauT_a.b = ", hyperParmTauT_a.b)
//------------------------------------------------------------------------------------------
/**
 * @object hyperParmTauT_b - contains parameters for the uniform prior of tauT_b
 * @property {number} a - value is the lower boundary of the uniform density
 * @property {number} b - value is the upper boundary of the uniform density
 */
var hyperParmTauT_b = {a:5.00, b:35.00}
console.log("hyperParmTauT_b.a = ", hyperParmTauT_b.a,
            ", hyperParmTauT_b.b = ", hyperParmTauT_b.b)
console.log("==============================================================================")
//------------------------------------------------------------------------------------------
// function definitions
//------------------------------------------------------------------------------------------
/**
 * @function runtime - method to compute the runtime in seconds and minutes
 */
var runTime = function() {
  var stopTime = Date.now()
  var runSecs = (stopTime - startTime)/1000
  var runMins = runSecs/60
  console.log("runtime in seconds = ", runSecs)
  console.log("runtime in minutes = ", runMins)}
//------------------------------------------------------------------------------------------
/**
 * @function oneSampleOfPrior - takes  o n e  sample from the uniform prior distributions
 * @returns {object} priorParamsObject - parameter object with properties
 *                                       priorTauT_a and priorTauT_b
 */
var oneSampleOfPrior = function () {
  var priorTauT_a = sample(Uniform({a:hyperParmTauT_a.a, b:hyperParmTauT_a.b}))
  var priorTauT_b = sample(Uniform({a:hyperParmTauT_b.a, b:hyperParmTauT_b.b}))
  var priorParamsObject = {priorTauT_a:priorTauT_a, priorTauT_b:priorTauT_b}
  return priorParamsObject
}
//------------------------------------------------------------------------------------------
/**
 * @description - Infer generates a bivariate uniform prior distribution for tauT
 * @variable {distribution} priorTauT - value is a WebPPL distribution object
 */
var priorTauT = Infer({model:oneSampleOfPrior, method: 'forward', samples: nTrials})
viz.marginals(priorTauT)
//------------------------------------------------------------------------------------------
/**
 * @description - marginalize generates a univariate uniform prior distribution for tauT_a
 * @variable {distribution} priorTauT_a - value is a WebPPL distribution object
 */
//------------------------------------------------------------------------------------------
var priorTauT_a = marginalize(priorTauT, 'priorTauT_a')
/**
 * @description - marginalize generates a univariate uniform prior distribution for tauT_b
 * @variable {distribution} priorTauT_b - value is a WebPPL distribution object
 */
var priorTauT_b = marginalize(priorTauT, 'priorTauT_b')
//------------------------------------------------------------------------------------------
/**
 * @function oneSampleOfModel - takes  o n e  sample from the posterior distribution
 * @returns {object} posteriorTauT_Object - posterior parameter object with properties
 *                                          'posteriorTauT_a' and 'posteriorTauT_b'
 */
var oneSampleOfModel = function() {
  /**
   * @variable {number} tauT_a - a sample from uniform prior distribution
   */
  var tauT_a = sample(priorTauT_a)
  /**
   * @variable {number} tauT_b - a sample from uniform prior distribution
   */
  var tauT_b = sample(priorTauT_b)
  //--------------------------------------------------------------------------------
  /**
   * @description - takes a datum as argument an returns a logLikelihood density
   *                of gamma distribution with shape and scale sampled from the prior
   * @function logLikelihood - value is logLikelihood_Gamma(datum;tauT_a, tauT_b)
   * @param {number} datum - one data value (in this case a RT)
   * @returns {number} logLikelihood_Gamma - logLikelihood_Gamma(datum;tauT_a, tauT_b)
   */
  var logLikelihood = function(datum) {
    var logLikelihood_Gamma = factor(Gamma({shape:tauT_a, scale:tauT_b}).score(datum))
    return logLikelihood_Gamma
  }
  //--------------------------------------------------------------------------------
  /**
  * @description sum(...) - value is logLikelihood function of data lnL(data; params)
  */
  sum(map(logLikelihood, data))   
  /**
  * @variable {object} posteriorTauT_Object - posterior parameter object with properties
  *                                          'posteriorTauT_a' and 'posteriorTauT_b'
  */
  var posteriorTauT_Object = {posteriorTauT_a:tauT_a, posteriorTauT_b:tauT_b}
  return posteriorTauT_Object}
//-----------------------------------------------------------------------------------
/**
 * @description - Infer generates the posterior distribution 'posteriorTauT'
 * @variable {distributionObject} posteriorTauT - bivariate posterior distribution
 */
var posteriorTauT = Infer({model:oneSampleOfModel, method:'MCMC', samples: nTrials,
                           burn:myBurnPeriod, lag:myLag})
viz.marginals(posteriorTauT)
//-----------------------------------------------------------------------------------
/**
 * @description - univariate marginal of bivariate posterior distribution of tauT
 * @variable {distributionObject} posteriorTauT_a - value is univariate posterior tauT_a
 */
var posteriorTauT_a = marginalize(posteriorTauT, 'posteriorTauT_a')
/**
 * @description - univariate marginal of bivariate posterior distribution of tauT
 * @variable {distributionObject} posteriorTauT_b - value is univariate posterior tauT_b
 */
var posteriorTauT_b = marginalize(posteriorTauT, 'posteriorTauT_b')
//-----------------------------------------------------------------------------------
/**
 * @description - Infer generates the 'posteriorPredictedTauT'-distribution
 *                this is the distribution of predicted response times taking parameter
 *                uncertainty into account
 * @function oneSampleOfPosteriorTauT - takes  o n e  sample from the posterior distributions
 *                                      posteriorTauT_a and posteriorTauT_b
 * @returns {object} posteriorTauT_Object - posterior parameter object with properties
 *                                          'shape' and 'scale'
 */
var oneSampleOfPosteriorTauT = function() {
  var sampledPosteriorTauT_a = sample(posteriorTauT_a)
  var sampledPosteriorTauT_b = sample(posteriorTauT_b)
  return sample(Gamma({shape:sampledPosteriorTauT_a, scale:sampledPosteriorTauT_b}))
}
var posteriorPredictedTauT = Infer({model:oneSampleOfPosteriorTauT, method: 'forward',
                                   samples:nTrials})
console.log("----------------------------------------------------------------------------")
console.log("posteriorPredictedTauT (this is the distribution of predicted response times")
console.log("======================  taking parameter uncertainty into account):")
viz.density(posteriorPredictedTauT)
/**
 * @description - descriptive statistics of a sample-generated gamma distribution
 * @function myTauXDistribution
 * @param {string} id - The identifier of the tauX distribution.
 * @param {distributionObject} tauXDistribution - tauX distribution (X = P, C, M, T)
 * @param {number} modeTauX - mode of tauX as a function of a and b
 *                            mode = (a-1)*b for a >= 1
 * @returns {object} meanSigmaTauObject - object with mean and sigma of TauX
 * @property {number} meanTauX - mean of tauX (X = P, C, M, T) or tau
 * @property {number} sigmaTauX - standard deviation of tauX (X = P, C, M, T) or tau
 */
var myTauXDescription = function(id, tauXDistribution, modeTauX) {
  var myTauXDistribution = { // extraction of probs and support from WebPPL tauX distribution
    probs: map(function(eventTuple){ // object to compute mean and sigma of tauX
      Math.exp(tauXDistribution.score(eventTuple))}, tauXDistribution.support()),
    support: tauXDistribution.support()}
  console.log("----------------------------------------------------------------------------")
  console.log(id)
  // mode(tauX), mean(tauX), variance(tauX) and sigma(tauX)
  console.log("mode = ", modeTauX)
  var meanTauX = sum(map2(function(value, prob) {
    value*prob},myTauXDistribution.support, myTauXDistribution.probs))
  console.log("mean = ", meanTauX)
  var sigmaTauX = Math.sqrt(sum(map2(function(value, prob) {
    Math.pow((value-meanTauX), 2)*prob},
                                     myTauXDistribution.support,
                                     myTauXDistribution.probs)))
  console.log("sigma = ", sigmaTauX)
  var tauX_Intval = {fast:meanTauX - nSigma * sigmaTauX, mean:meanTauX,
                     slow:meanTauX + nSigma * sigmaTauX}
  console.log("model-generated nSigma*sigma tau-interval: "); print(tauX_Intval)
  var meanSigmaTauXObject = {meanTauX: meanTauX, sigmaTauX: sigmaTauX}
  return meanSigmaTauXObject}
//
var xyz = myTauXDescription("posteriorPredictedTau", posteriorPredictedTauT, "unknown")
console.log("posteriorPredictedTauT, {bounds:[0,1500]}:")
viz.density(posteriorPredictedTauT, {bounds:[0,1500]})
console.log("posteriorPredictedTauT, {bounds:[0,1000]}:")
viz.density(posteriorPredictedTauT, {bounds:[0,1000]})
console.log("posteriorPredictedTauT, {bounds:[100, 500]}:")
viz.density(posteriorPredictedTauT, {bounds:[150, 500]})
//-----------------------------------------------------------------------------------
runTime()
console.log("==============================================================================")



WebPPL - Output

==============================================================================

PCM20180306_PersonalBayesianSimpleReactionTimeModel *** version 2018/03/27 ***

see also Simple Reaction Time, Example 9, Card, Moran & Newell, 1983, p.66    

see also https://www.humanbenchmark.com/tests/reactiontime/statistics 

==============================================================================

Input parameter:

nTrials =  20000

nSigma =  1

length of burn-in period =  2000

length of lag =  10

response time data = [ 458,292,228,403,271,420,350,235,260,306 ]

hyperParmTauT_a.a =  5 , hyperParmTauT_a.b =  50

hyperParmTauT_b.a =  5 , hyperParmTauT_b.b =  35

==============================================================================

 

priorTauT_a:

-------------------

Fig01 (see below)

-------------------

 

priorTauT_b:

-------------------

Fig02 (see below)

-------------------

 

posteriorTauT_a:

-------------------

Fig03 (see below)

-------------------

 

posteriorTauT_b:

-------------------

Fig04 (see below)

-------------------

----------------------------------------------------------------------------

posteriorPredictedTauT (this is the distribution of predicted response times

======================  taking parameter uncertainty into account):

-------------------

Fig05 (see below)

-------------------

----------------------------------------------------------------------------

posteriorPredictedTau

mode =  unknown

mean =  370.3064469388857

sigma =  218.18671171521027

model-generated nSigma*sigma tau-interval: 

{"fast":152.11973522367543,"mean":370.3064469388857,"slow":588.493158654096}

 

posteriorPredictedTauT, {bounds:[0,1500]}:

-------------------

Fig06 (see below)

-------------------

 

posteriorPredictedTauT, {bounds:[0,1000]}:

-------------------

Fig07 (see below)

-------------------

 

posteriorPredictedTauT, {bounds:[100, 500]}:

-------------------

Fig08 (see below)

-------------------

runtime in seconds =  1723.913

runtime in minutes =  28.731883333333332

==============================================================================
Weyfqbmaster (man/uuelc55a.wuesqraltefeld@uol.d5ue1jqn) (Changed: 2020-01-23)