prosafebeef Risk Model
Meat Safety Risk Assessment Model Tool
Page unquestionable under construction
This is supposed to be very interactive dashboard like demonstration of the different scenatio option that we can use for the Monte Carlo Risk simulation.
However, doing so in a SinglePageApplication with multiple javascript libraries turned out to be too much for single programmer...
Just to give you an impression of the type of work required here some code modules that are finished.
riskmodel.js
// javascript document with the model items and monte carlo simulation of the risk project // Samples object class which contains the actual simulation in time and during trasnformation and time for a given scenario // Scenario object holds the following subclasses : Pathogen, Production, Coldchain, Consumption, which in turn contain Dist(ributions) // Both the scenario it subclasses have a choise of predefined named setting but can also be custom configured during creation and later on. // erik.esveld@online.nl 2012 // version = 2 ///"use strict"; // samples object class creator takes a scenarion instance and number of samples // can perform calculations for a time or a specific moment function Samples(ascen, n) { // just store the references this.scenario = ascen; // number of samples this.n = n; // cooking time is just an fixed time for animation this.cookingtime = 10; // just because it is handy this.laglim = ascen.pathogen.lagwork / Math.LN10; // just because it is handy this.maxtime = ascen.shelflife * 24 + this.cookingtime; // inital loading this.logN0 = ascen.logN0.random(n); // transport time at 'Transported','Sold','Stored','Cooked' this.time = []; this.time[0] = ascen.transporttime.random(n) // make sure that the retail time is scaled such that its upper distribution limits = shelflife-transportime this.time[1] = ascen.retailtime.random(n); for (var i = 0; i < n; i++) this.time[1][i] *= (ascen.shelflife * 24 - this.time[0][i]) / ascen.retailtime.limits[1]; // storage time at home is random between shelflife - (retailtime + transporttime) this.time[2] = randoms(n); for (var i = 0; i < n; i++) this.time[2][i] *= ascen.shelflife * 24 - (this.time[0][i] + this.time[1][i]); // growth at 'Transported','Sold','Stored' no lag phase yet this.growth = []; this.growth[0] = ascen.transporttemp.random(n).map(ascen.pathogen.mumax, apat); this.growth[1] = ascen.retailtemp.random(n).map(ascen.pathogen.mumax, apat); this.growth[2] = ascen.storagetemp.random(n).map(ascen.pathogen.mumax, apat); // calculate growth from the mumax stored growth for (var j = 0; j < 3; j++) this.growth[j] = arraymultiply(this.growth[j], this.time[j]).map(mapscale(1 / Math.LN10)); // the growth during cooking is = -reductionrate/ln10*1 minute this.growth[3] = ascen.cookingtemp.random(n).map(ascen.pathogen.reductionrate, apat) .map(mapscale(-1 / Math.LN10)); // make grow times cumulative this.time[1] = arrayadd(this.time[0], this.time[1]); this.time[2] = arrayadd(this.time[1], this.time[2]); } // returns array of stages at time // transport starts at t=0 so time = -1 give 0=production // 'Production','Transport','Retail','Storage','Cooking','Consumption']; Samples.prototype.stage = function (t) { var p = []; for (var p = [], i = 0; i < this.n; i++) p.push(upperbin([-1, 0, this.time[0][i], this.time[1][i], this.time[2][i], this.time[2][i] + this.cookingtime], t)); return p; } Samples.prototype.attime = function (t) // return the contamination and stage at a time t // logN = logN0 + max(growth(t)-laglim,0) - reduction*(t-t_endofstorage)/t_cooking { var g, r; var logN = []; for (var i = 0; i < this.n; i++) { g = interp([0, this.time[0][i], this.time[1][i], this.time[2][i]], [0, this.growth[0][i], this.growth[0][i] + this.growth[1][i], this.growth[0][i] + this.growth[1][i] + this.growth[2][i]], t); r = Math.min(Math.max(t - this.time[2][i], 0), this.cookingtime) / this.cookingtime * this.growth[3][i]; logN.push(this.logN0[i] + Math.max(g - this.laglim, 0) + r); } var s = this.stage(t); var out = []; for (var i = 0; i < this.n; i++) out.push([logN[i], s[i]]); return out; } // return the contamination at a moment (e.g. 2 Sold) // 'Produced','Transported','Sold','Stored','Cooked' // logN = logN0 + max(growth-laglim,0) + reduction (which is negative) Samples.prototype.logN = function (moment) { switch (moment) { case 0: return this.logN0.slice(0); case 1: return arrayadd( this.logN0, this.growth[0] .map(mapshift(-this.laglim)) .map(mapmax(0)) ); case 2: return arrayadd( this.logN0, arrayadd( this.growth[0], this.growth[1] ) .map(mapshift(-this.laglim)) .map(mapmax(0)) ); case 3: return arrayadd( this.logN0, arrayadd( arrayadd( this.growth[0], this.growth[1] ), this.growth[2] ) .map(mapshift(-this.laglim)) .map(mapmax(0)) ); default: return arrayadd( this.growth[3], arrayadd( this.logN0, arrayadd( arrayadd( this.growth[0], this.growth[1] ), this.growth[2] ) .map(mapshift(-this.laglim)) .map(mapmax(0)) ) ); } } // return illness probabilities of the samples Samples.prototype.illnessprob = function () { var logNf = this.logNatmoment(4); var ip = []; for (var i = 0; i < this.n; i++) ip.push(this.pathogen.illnessprob(logNf[i])); return ip; } // quantile of samples (very slow) Samples.prototype.quantile = function (q) { return d3.quantile(this.sort(d3.ascending), q); } // scenarion object class contains all model settings // see Dist constructor for calling details function Scenario(ascen) { this.name = '-'; this.pathogen = new Pathogen(); this.production = new Production(); this.storage = new Storage(); this.consumption = new Consumption(); this.deflist = { }; this.typelist = []; // fill the typelist and add name deflist for (var lit in this.deflist) { this.typelist.push(lit); this.deflist[lit].name = "lit"; }; if (!!ascen) { // if agument was given if (!!this.deflist[ascen]) { // and agument is id of deflist ascen= this.deflist[ascen]; // then set to the default }; if (!!ascen.name) { this.name = ascen.name;} if (!!ascen.pathogen) { this.pathogen.constructor(ascen.pathogen);} if (!!ascen.production) { this.production.constructor(ascen.production);} if (!!ascen.storage) { this.production.constructor(ascen.production);} if (!!ascen.consumption) { this.consumption.constructor(ascen.consumption);} } return this; }; //Pathogen object class contains related properties and functions // see Dist creator for calling details function Pathogen(apat) { this.name = "-"; // typical concentration in feces which is potentially contaminating the carcas this.cfeces = new Dist({name: "Concentration in feces",unit: "log CFU/g",type: "Cumulative", points: [-1, 2, 3, 4, 5], cums: [0, 0.48, 0.54, 0.9, 1] }); // growth rate parameters // exponential growth rate maxmax = k*max(T-Tmin,0)^2 essentially same as equation Kostas // growth deltaLogN = max(mumax*t-h0,0)/LN10 (t since begin) this.k = 0.001; // hr-1 K-2 this.Tmin = 3; // degC this.h0= new Dist({type: "Normal", mu: 1.5, sigma: 0.5, limits: [0,3]}); this.mumax = function (T) { return this.k * Math.max(T - this.Tmin, 0) * Math.max(T - this.Tmin, 0); }; // exponential inactivation rate at temperature // deltaLogN = -reductionrate(T)/ln(10) *deltat in minutes // reductionrate(72) = ln(10)/D72 // reductionrate(T) = exp((T-72)/Z*ln(10))*reductionrate(72) this.D72 = 0.2; // minutes decimal reduction time at 72 this.Z = 10; // 10x decrease of decimal time every Z degrees higher than 72 this.reductionrate = function (T) {return Math.exp((T-72) / this.Z * Math.LN10) * Math.LN10 / this.D72;} // returns probability of illness after ingestion of logN (CFU/product) // CFD of normal distribution based on typical 5% and 50% quantiles fitting equation kostas // standard q50%-q5% = 1.65 so we apply a scaling to match with the given parameters this.illnessprob = function (logN) { return cnormal((logN-this.ID50)/(this.ID50 - this.ID5)*1.65);}; this.ID50 = 3; // 50% of consumers get ill after ingestion of ID50 logCFU this.ID5 = 2; // 5% of consumers get ill after ingestion of ID5 logCFU // default choises with only virulence parameters different // would be realy nice to include more kinds such as Salmonella, Listeria and C. Botulinum this.deflist = { EcoliO157: { name: "E-coli O157", ID50: 5.78, ID5: 4.33 }, EcoliO157Cassin: { name: "E-coli O157 Virulent", ID50: 3.4, ID5: 1.6 }, }; this.typelist = []; // fill the typelist and add name deflist for (var lit in this.deflist) { this.typelist.push(lit); }; if (!!apat) { // if agument was given if (!!this.deflist[apat]) { // and agument is id of deflist ascen= this.deflist[apat]; // then set to the default }; // copy the atrributes of the instance if given as first argument if (!!apat.name) { this.name = apat.name; } if (!!apat.cfeces) { this.cfeces.constructor(apat.cfeces); } if (!!apat.k) { this.k = apat.k; } if (!!apat.Tmin) { this.Tmin = apat.Tmin; } if (!!apat.h0) { this.h0.constructor(apat.h0); } if (!!apat.D72) { this.D72 = apat.D72; } if (!!apat.Z) { this.Z = apat.Z; } if (!!apat.ID50) { this.ID50 = apat.ID50; } if (!!apat.ID5) { this.ID5 = apat.ID5; } } return this; } // Production object encapsulates the related properties and function // see Dist creator for calling details function Production(aprod) { this.name = "-"; // determines the prevalence and error margin that a carcas is contaminated this.checked = 1; this.positive = 1; // determines if contaminated how much feces was spread on the carcas this.fecestransfer= new Dist({name: "feces transfer", unit: "log g/cm2", type: "Normal", mu : -5.1, sigma : 0.9, limits: [-9,-1]}); this.trimmingsurface = new Dist({name: "trimming surface", unit: "cm/g", type: "Uniform", limits: [0.25,1]}); // decontamination efficiency can be choosen this.decontamination = new Dist({name: 'none', unit: "-log", type: "Constant", mu: 0, limits : [-1,1]}); // default choises only country and decontamination this.deflist = { Ireland: { name : 'Ireland', checked : 332, positive: 7, decontamination: {name: 'none', unit: "-log", type: "Constant", mu: 0, limits : [-1,1]} }, Brazil: { name : 'Brazil', checked : 384, positive: 6, decontamination: {name: 'none', unit: "-log", type: "Constant", mu: 0, limits : [-1,1]} }, Greece: { name : 'Greece' , checked : 150, positive: 6, decontamination: {name: 'none', unit: "-log", type: "Constant", mu: 0, limits : [-1,1]} }, Poland: { name : 'Poland', checked : 486, positive: 80, decontamination: {name: 'none', unit: "-log", type: "Constant", mu: 0, limits : [-1,1]} }, PolandDecont: { name : 'Poland+Decontamination', checked : 486, positive: 80, decontamination: {name: 'decontamination', unit: "-log", type: "Normal", mu: 1.5, sigma : 0.5, limits : [0,3]} }, }; this.typelist = []; // fill the typelist and add name deflist for (var lit in this.deflist) { this.typelist.push(lit); }; if (!!aprod) { // if agument was given if (!!this.deflist[aprod]) { // and agument is id of deflist ascen= this.deflist[aprod]; // then set to the default }; // copy the relevant atrributes if defined in argument object if (!!aprod.name) { this.name = aprod.name; } if (!!aprod.checked) { this.checked = aprod.checked; } if (!!aprod.positive) { this.positive = aprod.positive; } if (!!aprod.fecestransfer) { this.fecestransfer.constructor(aprod.fecestransfer); } if (!!aprod.trimmingsurface) { this.trimmingsurface.constructor(aprod.trimmingsurface); } if (!!aprod.decontamination) { this.decontamination.constructor(aprod.decontamination); } } return this; }; // Cold chain object encapsulates the related properties and function // see Dist creator for calling details function Coldchain(acold) { this.name = "-"; // retail (package) size determines not so much except Poisson chance of having some bugs in it this.retailmass = new Dist({ name: 'Retail package', unit: 'gram', type: 'Triangular', limits: [300, 1000], peak: 500 }); // Shelf life is very important and detemines time in retail and domestic fridge this.shelflife = 7; //days // retail time (between production and sell ) needs to be multipled with shelflife this.retailtime = new Dist({ name: 'Retail time', unit: 'shelflife%', type: 'LogNormal', mu: 0.5 , sigma: 0.5 , limits: [1e-6, 1]}); // domestic time (between buy and eat) needs to be multipled with (shelflife - retailtime) this.hometime = new Dist({ name: 'Domestic time', unit: 'remainder%', type: 'Uniform', limits: [0, 1] }); // temperatures in retail and home this.retailtemp = new Dist({ name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }); this.hometemp = new Dist({ name: 'Home Temperature', unit: '°C', type: 'Normal', mu: 7.5, sigma: 3, limits: [0, 15]}); // default choises determine only shelf life and temperatures in retail // temperature higher than mu = 3.2 sigma = 2 are just demonstration estimations this.deflist = { Wholesale: { name: 'Wholesale', shelflife: 8, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }, }, Foodservice: { name: 'Food Service', shelflife: 7, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }, }, Supermarket: { name: 'Supermarket', shelflife: 6, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 4.2, sigma: 2.5, limits: [0, 7] }, }, Localshop: { name: 'Localshop', shelflife: 5, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 6, sigma: 3, limits: [0, 7] }, }, Market: { name: 'Market', shelflife: 3, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 8, sigma: 4, limits: [0, 7] }, }, }; this.typelist = []; // fill the typelist and add name deflist for (var lit in this.deflist) { this.typelist.push(lit); }; if (!!acold) { // if agument was given if (!!this.deflist[acold]) { // and agument is id of deflist ascen = this.deflist[acold]; // then set to the default }; // copy the atrributes of the instance if given as first argument if (!!acold.name) { this.name = acold.name; } if (!!acold.retailmass) { this.retailmass.constructor(acold.retailmass); } if (!!acold.shelflife) { this.shelflife = acold.shelflife; } if (!!acold.retailtime) { this.retailtime.constructor(acold.retailtime); } if (!!acold.hometime) { this.hometime.constructor(acold.hometime); } if (!!acold.retailtemp) { this.retailtemp.constructor(acold.retailtemp); } if (!!acold.hometemp) { this.hometemp.constructor(acold.hometemp); } } return this; }; // Cold chain object encapsulates the related properties and function // see Dist creator for calling details function Coldchain(acold) { this.name = "-"; // retail (package) size determines not so much except Poisson chance of having some bugs in it this.retailmass = new Dist({ name: 'Retail package', unit: 'gram', type: 'Triangular', limits: [300, 1000], peak: 500 }); // Shelf life is very important and detemines time in retail and domestic fridge this.shelflife = 7; //days // retail time (between production and sell ) needs to be multipled with shelflife this.retailtime = new Dist({ name: 'Retail time', unit: 'shelflife%', type: 'LogNormal', mu: 0.5, sigma: 0.5, limits: [1e-6, 1] }); // domestic time (between buy and eat) needs to be multipled with (shelflife - retailtime) this.hometime = new Dist({ name: 'Domestic time', unit: 'remainder%', type: 'Uniform', limits: [0, 1] }); // temperatures in retail and home this.retailtemp = new Dist({ name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }); this.hometemp = new Dist({ name: 'Home Temperature', unit: '°C', type: 'Normal', mu: 7.5, sigma: 3, limits: [0, 15] }); // default choises determine only shelf life and temperatures in retail // temperature higher than mu = 3.2 sigma = 2 are just demonstration estimations this.deflist = { Wholesale: { name: 'Wholesale', shelflife: 8, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }, }, Foodservice: { name: 'Food Service', shelflife: 7, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 3.2, sigma: 2, limits: [0, 7] }, }, Supermarket: { name: 'Supermarket', shelflife: 6, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 4.2, sigma: 2.5, limits: [0, 7] }, }, Localshop: { name: 'Localshop', shelflife: 5, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 6, sigma: 3, limits: [0, 7] }, }, Market: { name: 'Market', shelflife: 3, retailtemp: { name: 'Retail Temperature', unit: '°C', type: 'Normal', mu: 8, sigma: 4, limits: [0, 7] }, }, }; this.typelist = []; // fill the typelist and add name deflist for (var lit in this.deflist) { this.typelist.push(lit); }; if (!!acold) { // if agument was given if (!!this.deflist[acold]) { // and agument is id of deflist ascen = this.deflist[acold]; // then set to the default }; // copy the atrributes of the instance if given as first argument if (!!acold.name) { this.name = acold.name; } if (!!acold.retailmass) { this.retailmass.constructor(acold.retailmass); } if (!!acold.shelflife) { this.shelflife = acold.shelflife; } if (!!acold.retailtime) { this.retailtime.constructor(acold.retailtime); } if (!!acold.hometime) { this.hometime.constructor(acold.hometime); } if (!!acold.retailtemp) { this.retailtemp.constructor(acold.retailtemp); } if (!!acold.hometemp) { this.hometemp.constructor(acold.hometemp); } } return this; };
riskdist.js
// javascript document with the statistical distributions of the risk project // defines Risk object class and calculation functions // erik.esveld@online.nl 2012 // version = 2 "use strict"; // Dist object class to deal with continues distributions // creator takes optionally as setting // an instance -> var adist = new Dist(); var bdist = new Dist(adist); // or object with literal parameters -> var adist = new Dist({type:'Uniform', limits:[0,2]}) // or string with a default id -> var adist = new Dist('Normal') // used to create and change the distributions // and to create random samples distributed accordingly function Dist(adist) { // properties this.name = "-"; this.unit = "-"; this.type = "Uniform"; this.limits = [0, 2]; // lognormal --> limits[0] >=0 Triangular -> edges this.mu = 1; //normal -> mean and lognormal -> 'location'or mean of ln(X) ' this.sigma = 0.5; // normal -> stdev and lognormal -> 'scale' or stdev of ln(X) this.peak = 1; // triangular this.points = [0, 1, 2]; // cumulative points acending ! this.cums = [0, 0.5, 1]; // corresponding cumulative chance inbetween 0-1 acending ! this.deflist = { Constant: { mu: 1, limits: [0, 2] }, Uniform: { limits: [0, 2] }, Normal: { limits: [0, 2], mu: 1, sigma: .5 }, LogNormal: { limits: [0, 5], mu: 1, sigma: 1 }, Triangular: { limits: [0, 2], peak: 1 }, Cumulative: { limits: [0, 2], points: [0, 1, 2], cums: [0, 0.5, 1] }, }; this.typelist = []; // fill the typelist and add name and type to deflist for (var lit in this.deflist) { this.typelist.push(lit); this.deflist[lit].name = "default"; this.deflist[lit].type = lit; }; // set the properties from adist // we do it like this because if we poll the item of adist we also get non wanted properties suchs the class methods if (!!adist) { // if argument is given if (!!this.deflist[adist]) { // and argument is id of deflist adist = this.deflist[adist]; // set to argument default } if (!!adist.name) { this.name = adist.name; } if (!!adist.unit) { this.unit = adist.unit; } if (!!adist.type) { this.type = adist.type; } if (!!adist.limits) { this.limits = adist.limits.slice(0); } if (!!adist.mu) { this.mu = adist.mu; } if (!!adist.sigma) { this.sigma = adist.sigma; } if (!!adist.peak) { this.peak = adist.peak; } if (!!adist.points) { this.points = adist.points; } // use the .creator if it was an object instead of simple array ; if (!!adist.cums) { this.cums = adist.cums; } } return this.paramcheck(); } // set type do some sanity checking and return this Dist.prototype.paramcheck = function () { this.limits.sort(numascending); switch (this.type) { case "Constant": this.limits = [this.mu - 1, this.mu + 1]; break; case "LogNormal": if (this.limits[0] < 0) { this.limits[0] *= -1 }; if (this.limits[1] < this.limits[0]) { this.limits[1] *= -1 }; break; case "Triangular": if (upperbin(this.limits, this.peak) != 1) { this.peak = 0.5 * (this.limits[0] + this.limits[1]); } break; case "Cumulative": this.cums[0] = 0; this.cums[length(this.points) - 1] = 1; this.points.sort(numascending); this.cums.sort(numascending); default: return this; } return this; } // Dist probability density function Dist.prototype.pdf = function (x) { if (x < this.limits[0] || x > this.limits[1]) { return 0; } this.paramcheck(); switch (this.type) { case "Constant": return (x == this.mu ? 1 : 0); case "Uniform": return 1 / (this.limits[1] - this.limits[0]); case "Normal": var cll = cnormal((this.limits[0] - this.mu) / this.sigma); var cul = cnormal((this.limits[1] - this.mu) / this.sigma); return this.pnormal((x - this.mu) / this.sigma) / this.sigma / (cul - cll); case "LogNormal": if (x <= 1e-6) { return 0; } else { var cll = cnormal(Math.log(Math.max((this.limits[0] - this.mu) / this.sigma, 1e-6))); var cul = cnormal(Math.log(Math.max((this.limits[1] - this.mu) / this.sigma), 2e-6)); return this.pnormal((Math.log(x) - this.mu) / this.sigma) / this.sigma / (cul - cll); } case "Triangular": return ptriangular(x, this.limits[0], this.limits[1], this.peak); case "Cumulative": return pcumulative(x, this.points, this.cums); default: return 1 / (this.limits[1] - this.limits[0]); } } // Dist cumulative density function Dist.prototype.cdf = function (x) { if (x < this.limits[0]) { return 0; } if (x > this.limits[1]) { return 1; } switch (this.type) { case "Constant": return (x < this.mu ? 0 : 1); case "Uniform": return (x - this.limits[0]) / (this.limits[1] - this.limits[0]); case "Normal": var cll = cnormal((this.limits[0] - this.mu) / this.sigma); var cul = cnormal((this.limits[1] - this.mu) / this.sigma); return (cnormal((x - this.mu) / this.sigma) - cll) / (cul - cll); case "LogNormal": if (x <= 1e-6) { return 0; } else { var cll = cnormal(Math.log(Math.max((this.limits[0] - this.mu) / this.sigma, 1e-6))); var cul = cnormal(Math.log(Math.max((this.limits[1] - this.mu) / this.sigma), 2e-6)); return (cnormal((Math.log(x) - this.mu) / this.sigma) - cll) / (cul - cll); } case "Triangular": return ctriangular(x, this.limits[0], this.limits[1], this.peak); case "Cumulative": return ccumulative(x, this.points, this.cums); default: return (x - this.limits[0]) / (this.limits[1] - this.limits[0]); } } // Dist inverse cumulative density function or quantile Dist.prototype.icdf = function (y) { if (y <= 0) { return this.limits[0] } if (y >= 1) { return this.limits[1] } switch (this.type) { case "Constant": return this.mu; case "Uniform": return this.limits[0] + y * (this.limits[1] - this.limits[0]); case "Normal": var cll = cnormal((this.limits[0] - this.mu) / this.sigma); var cul = cnormal((this.limits[1] - this.mu) / this.sigma); return icnormal(y * (cul - cll) + cll) * this.sigma + this.mu; case "LogNormal": if (x <= 1e-6) { return 0; } else { var cll = cnormal(Math.log(Math.max((this.limits[0] - this.mu) / this.sigma, 1e-6))); var cul = cnormal(Math.log(Math.max((this.limits[1] - this.mu) / this.sigma), 2e-6)); return Math.exp(icnormal(y * (cul - cll) + cll) * this.sigma + this.mu); } case "Triangular": return ictriangular(y, this.limits[0], this.limits[1], this.peak); case "Cumulative": return iccumulative(y, this.points, this.cums); default: return (x - this.limits[0]) / (this.limits[1] - this.limits[0]); } } // returns array of random values distributed according to the distribution type // output lies between limits Dist.prototype.random = function (n) { switch (this.type) { case "Constant": return ones(n).map(mapscale(this.mu)); case "Uniform": return randoms(n).map(maprange(this.limits[0], this.limits[1])); case "Normal": var z = zeros(n); for (var i = 0; i < n; i++) { do { z[i] = rnormal() * this.sigma + this.mu; } while (z[i] < this.limits[0] || z[i] > this.limits[1]) } return z; case "LogNormal": var z = zeros(n); for (var i = 0; i < length(z) ; i++) { do { z[i] = Math.exp(rnormal() * this.sigma + this.mu); } while (z[i] < this.limits[0] || z[i] > this.limits[1]) } return z; case "Triangular": var z = random(n); for (var i = 0; i < n; i++) { z[i] = ictrangular(z[i], this.limits[0], this.limits[1], this.peak); } return z; case "Cumulative": var z = random(n); for (var i = 0; i < n; i++) { z[i] = icumulative(z[i], this.points, this.cums); } return z; default: return randoms(n).map(maprange(this.limits[0], this.limits[1])); } } /////////// DISTRIBUTION FUNCTIONS /////////// // normal distribution functions //PDF function pnormal(x) { return 0.4 * Math.exp(-x * x / 2); } //CDF function cnormal(x) { return 0.5 * (1 + erf(x / Math.SQRT2)); } //ICDF (0<=y<=1) function icnormal(y) { return Math.SQRT2 * ierf(2 * y - 1); } // returns a normal random values Box Muller method function rnormal() { var r = Math.ln(-2 * Math.random()); var alpha = 2 * Math.PI * Math.random(); return r * Math.cos(alpha); } // triangular distrubution functions //PDF function ptriangular(x, ll, ul, mx) { if (x < ll || x > ul) { return 0; } if (x < mx) { return 2 * (x - ll) / (ul - ll) / (mx - ll); } else { return 2 * (ul - x) / (ul - ll) / (ul - mx); } } // CDF function ctriangular(x, ll, ul, mx) { if (x < ll) { return 0; } if (x > ul) { return 1; } if (x < mx) { return (x - ll) * (x - ll) / (ul - ll) / (mx - ll); } else { return 1 - ((ul - x) * (ul - x) / (ul - ll) * (ul - mx)); } } // ICDF (0<=y<=1) function ictrangular(y, ll, ul, mx) { if (y < ((mx - ll) / (ul - ll))) { return ll + Math.sqrt(y * (mx - ll) * (ul - ll)); } else { return ul - Math.sqrt((1 - y) * (ul - mx) * (ul - ll)); } } // cumulatice piecewise distributionfunctions // PDF function pcumulative(x, points, cums) { var n = length(points); var p = zeros[n + 1]; for (var i = 0; i < n - 1; i++) { p[i + 1] = (cums[i + 1] - cums[i]) / (points[i + 1] - points[i]); } return p[upperbin(points, x)]; } // CDF function ccumulative(x, points, cums) { return interp(points, cums, x); } // ICDF function iccumulative(y, points, cums) { return interp(cums, points, y); } // discrete Poisson functions // Create a Poison ditributed random variable function rpoisson(lambda) { if (lambda >= 30.0) { return poissonLargeAprox(lambda); // return poissonLarge(lambda) } else { // Algorithm due to Donald Knuth, 1969. var L = exp(-lambda); var p = Math.random(); var k = 0; while (p > L) { p *= Math.random(); k++; } return k; } } // Create a Poison ditributed random variable Greater Than 0 function rpoissonGT0(lambda) { var k = 0; if (lambda >= 30.0) { do { k = poissonLargeAprox(lambda); // var k = poissonLarge(lambda) } while (k == 0); return k; } else { // Algorithm due to Donald Knuth, 1969. var L = exp(-lambda); // p > L for sure var p = 1 - Math.random() * (1 - L); while (p > L) { p *= Math.random(); k++; } return k; } } // Probility for a number Greater Than 0 is 1-exp(-lambda) function ppoissonGT0(lambda) { return 1 - Math.exp(-lambda); } // crude but fast approximation based on normal distribution and continuety correction function poissonLargeAprox(lambda) { return Math.abs(Math.floor(rnormal() * Math.sqrt(lambda) + lambda + 0.5)); } // from http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values // "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson // Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979) // The article is on pages 29-35. The algorithm given here is on page 32. function poissonLarge(lambda) { var c = 0.767 - 3.36 / lambda; var beta = Math.PI / Math.sqrt(3.0 * lambda); var alpha = beta * lambda; var u, x, n, v, y, temp, lhs, rhs; var k = Math.log(c) - lambda - Math.log(beta); do { do { u = Math.random(); x = (alpha - Math.log((1.0 - u) / u)) / beta; n = Math.floor(x + 0.5); } while (n < 0); v = Math.random(); y = alpha - beta * x; temp = 1.0 + Math.exp(y); lhs = y + Math.log(v / (temp * temp)); rhs = k + n * Math.log(lambda) - LogFactorial(n); } while (lhs <= rhs); return n; } /////////////// UTILITY FUCTIONS //////////////////// // find first index for which B[i]>=a // returns 0 if B[0]>=a // returns B.length as B[last]a function upperbin(B, a) { var i = 0; while (B[i] < a) { ++i; } return i; } // interpolation of a value on linear piecewise function // flat extrapolation outside X domain function interp(X, Y, xi) { // X and Y must be same length arrays // X is monotonic increasing // xi is scalar with min(X) <= xi <= max(X) // first find the index for which X[i-1] < xi < X[i] var i = upperbin(X, xi); if (i == 0) return Y[0]; if (i == X.length) return Y[i - 1]; if ((X[i] - X[i - 1]) != 0) //maybe we shoudl leave this out return Y[i - 1] + (Y[i] - Y[i - 1]) * (xi - X[i - 1]) / (X[i] - X[i - 1]); else return (Y[i - 1] + Y[i]) / 2; } // integration of stepwise function starting at 0 // int= 0 for xi <0 // int = int(max(X)) for xi > max(X) // is not the same as cdf since freq[i] is already the intergral of the probability over bin[i-1] to bin[i] function stepintegral(X, Y, xi) { // X and Y must be same length arrays // X is monotonic increasing // xi is scalar with min(X) <= xi <= max(X) // first find the index for which X[i-1] < xi < X[i] if (xi <= 0) return 0; var j = upperbin(X, xi); if (j == X.length) { j = X.length - 1; xi = X[j]; } var int = 0; // integration of the whole sections for (var i = 0; i < j; i++) { int += Y[i] * (X[i] - (X[i - 1] || 0)); } // add the integral of the partial last section int += Y[j] * (xi - (X[j - 1] || 0)); return int; } // integration of linear piecewise function // int= 0 for ximax(X) function linintegral(X, Y, xi) { // X and Y must be same length arrays // X is monotonic increasing // xi is scalar with min(X) <= xi <= max(X) // first find the index for which X[i-1] < xi < X[i] var j = upperbin(X, xi); if (j == 0) return 0; if (j == X.length) { j = X.length - 1; xi = X[j]; } var int = 0; // integration of the whole sections for (var i = 0; i < (j - 1) ; i++) { int += (Y[i] + Y[i + 1]) * (X[i + 1] - X[i]) / 2; } // add the integral of the partial last section var yi = interp([X[j - 1], X[j]], [Y[j - 1], Y[j]], xi); int += (Y[j - 1] + yi) * (xi - X[j - 1]) / 2; return int; } // array creation functions // returns an array of length n with zeros function zeros(n) { var a = []; for (var i = 0; i < n; i++) a.push(0); return a; } // returns an array of length n with ones function ones(n) { var a = []; for (var i = 0; i < n; i++) a.push(1); return a; } // returns an array of length n with random values between 0 and <1 function randoms(n) { var a = []; for (var i = 0; i < n; i++) a.push(Math.random()); return a; } // array functions // return boolean if araya contains a function arraycontains(arraya, a) { return (arraya.indexOf(a) != -1); } // return sum arraya function arraysum(arraya) { var a = 0; for (var i = 0; i < arraya.length; i++) a += arraya[i]; return a; } // returns cumulativce sum arraya function arraycumsum(arraya) { var a = arraya[0]; for (var i = 1; i < arraya.length; i++) a.push(a[i - 1] + arraya[i]); return a; } // returns arraya .+ arrayb function arrayadd(arraya, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(arraya[i] + arrayb[i]); return a; } // returns arraya .- arrayb function arrayminus(arraya, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(arraya[i] - arrayb[i]); return a; } // returns arraya .* arrayb function arraymultiply(arraya, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(arraya[i] * arrayb[i]); return a; } // returns arraya ./ arrayb NO CHECK for division by zero function arraydivision(arraya, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(arraya[i] / arrayb[i]); return a; } // returns min(arraya,arrayb) function arraymin(array, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(Math.min(arraya[i], arrayb[i])); return a; } // returns max(arraya,arrayb) function arraymax(array, arrayb) { var a = []; for (var i = 0; i < arraya.length; i++) a.push(Math.max(arraya[i], arrayb[i])); return a; } // ascessor functions for array.map() function mapscale(s) { return function (d) { return d * s } } function mapshift(s) { return function (d) { return d + s } } function maprange(l, u) { return function (d) { return d * (u - l) + l } } function mapinterp(X, Y) { return function (d) { return interp(X, Y, d) } } function mapkey(key) { return function (d) { return d[key] } } function mapmin(s) { return function (d) { return Math.min(d, s) } } function mapmax(s) { return function (d) { return Math.max(d, s) } } // arraysort function function numascending(a, b) { return a - b }; ////////////////////////// Extra math functions ///////////////////////////////////// function sign(x) { return x > 0 ? 1 : x == 0 ? 0 : -1; } // error function approximation function erf(x) { return sign(x) * Math.sqrt(1 - Math.exp(-x * x * (4 / Math.PI + 0.14 * x * x) / (1 + 0.14 * x * x))); } // inverse error function approximation function ierf(y) { var a = Math.ln(1 - y * y); var b = 2 / 0.14 / Math.PI + a / 2; return sign(y) * Math.sqrt(Math.sqrt(b * b - a / 0.14) - b); } function factorial(n) { return (n < 2) ? 1 : n * factorial(n - 1); } // from http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values function LogFactorial(n) { if (n > 254) { var x = n + 1; return (x - 0.5) * Math.log(x) - x + 0.5 * Math.log(2 * Math.PI) + 1.0 / (12.0 * x); } else { var lf = [ 0.000000000000000, 0.000000000000000, 0.693147180559945, 1.791759469228055, 3.178053830347946, 4.787491742782046, 6.579251212010101, 8.525161361065415, 10.604602902745251, 12.801827480081469, 15.104412573075516, 17.502307845873887, 19.987214495661885, 22.552163853123421, 25.191221182738683, 27.899271383840894, 30.671860106080675, 33.505073450136891, 36.395445208033053, 39.339884187199495, 42.335616460753485, 45.380138898476908, 48.471181351835227, 51.606675567764377, 54.784729398112319, 58.003605222980518, 61.261701761002001, 64.557538627006323, 67.889743137181526, 71.257038967168000, 74.658236348830158, 78.092223553315307, 81.557959456115029, 85.054467017581516, 88.580827542197682, 92.136175603687079, 95.719694542143202, 99.330612454787428, 102.968198614513810, 106.631760260643450, 110.320639714757390, 114.034211781461690, 117.771881399745060, 121.533081515438640, 125.317271149356880, 129.123933639127240, 132.952575035616290, 136.802722637326350, 140.673923648234250, 144.565743946344900, 148.477766951773020, 152.409592584497350, 156.360836303078800, 160.331128216630930, 164.320112263195170, 168.327445448427650, 172.352797139162820, 176.395848406997370, 180.456291417543780, 184.533828861449510, 188.628173423671600, 192.739047287844900, 196.866181672889980, 201.009316399281570, 205.168199482641200, 209.342586752536820, 213.532241494563270, 217.736934113954250, 221.956441819130360, 226.190548323727570, 230.439043565776930, 234.701723442818260, 238.978389561834350, 243.268849002982730, 247.572914096186910, 251.890402209723190, 256.221135550009480, 260.564940971863220, 264.921649798552780, 269.291097651019810, 273.673124285693690, 278.067573440366120, 282.474292687630400, 286.893133295426990, 291.323950094270290, 295.766601350760600, 300.220948647014100, 304.686856765668720, 309.164193580146900, 313.652829949878990, 318.152639620209300, 322.663499126726210, 327.185287703775200, 331.717887196928470, 336.261181979198450, 340.815058870798960, 345.379407062266860, 349.954118040770250, 354.539085519440790, 359.134205369575340, 363.739375555563470, 368.354496072404690, 372.979468885689020, 377.614197873918670, 382.258588773060010, 386.912549123217560, 391.575988217329610, 396.248817051791490, 400.930948278915760, 405.622296161144900, 410.322776526937280, 415.032306728249580, 419.750805599544780, 424.478193418257090, 429.214391866651570, 433.959323995014870, 438.712914186121170, 443.475088120918940, 448.245772745384610, 453.024896238496130, 457.812387981278110, 462.608178526874890, 467.412199571608080, 472.224383926980520, 477.044665492585580, 481.872979229887900, 486.709261136839360, 491.553448223298010, 496.405478487217580, 501.265290891579240, 506.132825342034830, 511.008022665236070, 515.890824587822520, 520.781173716044240, 525.679013515995050, 530.584288294433580, 535.496943180169520, 540.416924105997740, 545.344177791154950, 550.278651724285620, 555.220294146894960, 560.169054037273100, 565.124881094874350, 570.087725725134190, 575.057539024710200, 580.034272767130800, 585.017879388839220, 590.008311975617860, 595.005524249382010, 600.009470555327430, 605.020105849423770, 610.037385686238740, 615.061266207084940, 620.091704128477430, 625.128656730891070, 630.172081847810200, 635.221937855059760, 640.278183660408100, 645.340778693435030, 650.409682895655240, 655.484856710889060, 660.566261075873510, 665.653857411105950, 670.747607611912710, 675.847474039736880, 680.953419513637530, 686.065407301994010, 691.183401114410800, 696.307365093814040, 701.437263808737160, 706.573062245787470, 711.714725802289990, 716.862220279103440, 722.015511873601330, 727.174567172815840, 732.339353146739310, 737.509837141777440, 742.685986874351220, 747.867770424643370, 753.055156230484160, 758.248113081374300, 763.446610112640200, 768.650616799717000, 773.860102952558460, 779.075038710167410, 784.295394535245690, 789.521141208958970, 794.752249825813460, 799.988691788643450, 805.230438803703120, 810.477462875863580, 815.729736303910160, 820.987231675937890, 826.249921864842800, 831.517780023906310, 836.790779582469900, 842.068894241700490, 847.352097970438420, 852.640365001133090, 857.933669825857460, 863.231987192405430, 868.535292100464630, 873.843559797865740, 879.156765776907600, 884.474885770751830, 889.797895749890240, 895.125771918679900, 900.458490711945270, 905.796028791646340, 911.138363043611210, 916.485470574328820, 921.837328707804890, 927.193914982476710, 932.555207148186240, 937.921183163208070, 943.291821191335660, 948.667099599019820, 954.046996952560450, 959.431492015349480, 964.820563745165940, 970.214191291518320, 975.612353993036210, 981.015031374908400, 986.422203146368590, 991.833849198223450, 997.249949600427840, 1002.670484599700300, 1008.095434617181700, 1013.524780246136200, 1018.958502249690200, 1024.396581558613400, 1029.838999269135500, 1035.285736640801600, 1040.736775094367400, 1046.192096209724900, 1051.651681723869200, 1057.115513528895000, 1062.583573670030100, 1068.055844343701400, 1073.532307895632800, 1079.012946818975000, 1084.497743752465600, 1089.986681478622400, 1095.479742921962700, 1100.976911147256000, 1106.478169357800900, 1111.983500893733000, 1117.492889230361000, 1123.006317976526100, 1128.523770872990800, 1134.045231790853000, 1139.570684729984800, 1145.100113817496100, 1150.633503306223700, 1156.170837573242400 ]; return lf[n]; } }