time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

beta-test.js (2478B)


      1 var vows = require('vows');
      2 var assert = require('assert');
      3 var suite = vows.describe('jStat.distribution');
      4 
      5 require('../env.js');
      6 
      7 suite.addBatch({
      8   'beta pdf': {
      9     'topic': function() {
     10       return jStat;
     11     },
     12     'check pdf calculation': function(jStat) {
     13       // Non-log form of the Beta pdf
     14       function pdf(x, alpha, beta) {
     15         if (x > 1 || x < 0)
     16           return 0;
     17         return (Math.pow(x, alpha - 1) * Math.pow(1 - x, beta - 1)) /
     18             jStat.betafn(alpha, beta);
     19       }
     20 
     21       var tol = 0.0000001;
     22       var args = [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1];
     23       var arg;
     24 
     25       for (var i = 0; i < args.length; i++) {
     26         arg = args[i];
     27         assert.epsilon(tol, jStat.beta.pdf(arg, 0.1, 0.1), pdf(arg, 0.1, 0.1));
     28         assert.epsilon(tol, jStat.beta.pdf(arg, 1, 1), pdf(arg, 1, 1));
     29         assert.epsilon(tol, jStat.beta.pdf(arg, 10, 50), pdf(arg, 10, 50));
     30 
     31         // Show that the log form of the pdf performs better for
     32         // large parameter values.
     33         assert(!isNaN(jStat.beta.pdf(arg, 1000, 5000)),
     34                'New Beta pdf is NaN for large parameter values.');
     35         assert(isNaN(pdf(arg, 1000, 5000)),
     36                'Old Beta pdf is not NaN for large parameter values.');
     37       }
     38 
     39       assert.equal(jStat.beta.pdf(0, 1, 4), 4);
     40       assert.equal(jStat.beta.pdf(1, 4, 1), 4);
     41       assert.equal(jStat.beta.pdf(0.5, 200, 4000), 0);
     42     },
     43     // checked against R code:
     44     //   options(digits=10)
     45     //   # Using mode definition from: https://en.wikipedia.org/wiki/Beta_distribution
     46     //   beta.mode <- function (a, b) {(a-1)/(a+b-2)}
     47     //   beta.mode(2.05, 2)
     48     //   beta.mode(5, 10)
     49     //   beta.mode(3, 3)
     50     'check mode calculation': function(jStat) {
     51       var tol = 0.0000001;
     52       assert.epsilon(tol, jStat.beta.mode(5, 10), 0.3076923077);
     53       assert.epsilon(tol, jStat.beta.mode(2.05, 2), 0.512195122);
     54       assert.epsilon(tol, jStat.beta.mode(3, 3), 0.5);
     55     },
     56     // checked against R's qbeta(p, shape1, shape2, ncp=0, lower.tail=TRUE, log.p=FALSE) from package 'stats':
     57     //   options(digits=10)
     58     //   qbeta(0.5, 5, 10)
     59     //   qbeta(0.5, 2.05, 2)
     60     //   qbeta(0.5, 3, 3)
     61     'check median calculation': function(jStat) {
     62       var tol = 0.0000001;
     63       assert.epsilon(tol, jStat.beta.median(5, 10), 0.3257511553);
     64       assert.epsilon(tol, jStat.beta.median(2.05, 2), 0.5072797399);
     65       assert.epsilon(tol, jStat.beta.median(3, 3), 0.5);
     66     }
     67   }
     68 });
     69 
     70 suite.export(module);