/* global module *//* global module */A simple, literate statistics system. The code below uses the
+Javascript module pattern,
+eventually assigning simple-statistics to ss in browsers or the
+exports object for node.js
(function() {
+ var ss = {};
+
+ if (typeof module !== 'undefined') { module.exports = ss;
+ } else {Otherwise, in a browser, we assign ss to the window object,
+so you can simply refer to it as ss.
this.ss = ss;
+ }Simple linear regression +is a simple way to find a fitted line +between a set of coordinates.
+ + function linear_regression() {
+ var linreg = {},
+ data = [];Assign data to the model. Data is assumed to be an array.
+ + linreg.data = function(x) {
+ if (!arguments.length) return data;
+ data = x.slice();
+ return linreg;
+ };Calculate the slope and y-intercept of the regression line +by calculating the least sum of squares
+ + linreg.mb = function() {
+ var m, b;Store data length in a local variable to reduce +repeated object property lookups
+ + var data_length = data.length;if there’s only one point, arbitrarily choose a slope of 0 +and a y-intercept of whatever the y of the initial point is
+ + if (data_length === 1) {
+ m = 0;
+ b = data[0][1];
+ } else {Initialize our sums and scope the m and b
+variables that define the line.
var sum_x = 0, sum_y = 0,
+ sum_xx = 0, sum_xy = 0;Use local variables to grab point values +with minimal object property lookups
+ + var point, x, y;Gather the sum of all x values, the sum of all +y values, and the sum of x^2 and (x*y) for each +value.
+In math notation, these would be SS_x, SS_y, SS_xx, and SS_xy
+ + for (var i = 0; i < data_length; i++) {
+ point = data[i];
+ x = point[0];
+ y = point[1];
+
+ sum_x += x;
+ sum_y += y;
+
+ sum_xx += x * x;
+ sum_xy += x * y;
+ }m is the slope of the regression line
m = ((data_length * sum_xy) - (sum_x * sum_y)) /
+ ((data_length * sum_xx) - (sum_x * sum_x));b is the y-intercept of the line.
b = (sum_y / data_length) - ((m * sum_x) / data_length);
+ }Return both values as an object.
+ + return { m: m, b: b };
+ };a shortcut for simply getting the slope of the regression line
+ + linreg.m = function() {
+ return linreg.mb().m;
+ };a shortcut for simply getting the y-intercept of the regression +line.
+ + linreg.b = function() {
+ return linreg.mb().b;
+ };This is called after .data() and returns the
+equation y = f(x) which gives the position
+of the regression line at each point in x.
linreg.line = function() {Get the slope, m, and y-intercept, b, of the line.
var mb = linreg.mb(),
+ m = mb.m,
+ b = mb.b;Return a function that computes a y value for each
+x value it is given, based on the values of b and a
+that we just computed.
return function(x) {
+ return b + (m * x);
+ };
+ };
+
+ return linreg;
+ }The r-squared value of data compared with a function f
+is the sum of the squared differences between the prediction
+and the actual value.
function r_squared(data, f) {
+ if (data.length < 2) return 1;Compute the average y value for the actual +data set in order to compute the +total sum of squares
+ + var sum = 0, average;
+ for (var i = 0; i < data.length; i++) {
+ sum += data[i][1];
+ }
+ average = sum / data.length;Compute the total sum of squares - the +squared difference between each point +and the average of all points.
+ + var sum_of_squares = 0;
+ for (var j = 0; j < data.length; j++) {
+ sum_of_squares += Math.pow(average - data[j][1], 2);
+ }Finally estimate the error: the squared +difference between the estimate and the actual data +value at each point.
+ + var err = 0;
+ for (var k = 0; k < data.length; k++) {
+ err += Math.pow(data[k][1] - f(data[k][0]), 2);
+ }As the error grows larger, its ratio to the +sum of squares increases and the r squared +value grows lower.
+ + return 1 - (err / sum_of_squares);
+ }This is a naïve bayesian classifier that takes +singly-nested objects.
+ + function bayesian() {The bayes_model object is what will be exposed
+by this closure, with all of its extended methods, and will
+have access to all scope variables, like total_count.
var bayes_model = {},The number of items that are currently +classified in the model
+ + total_count = 0,Every item classified in the model
+ + data = {};Train the classifier with a new item, which has a single +dimension of Javascript literal keys and values.
+ + bayes_model.train = function(item, category) {If the data object doesn’t have any values +for this category, create a new object for it.
+ + if (!data[category]) data[category] = {};Iterate through each key in the item.
+ + for (var k in item) {
+ var v = item[k];Initialize the nested object data[category][k][item[k]]
+with an object of keys that equal 0.
if (data[category][k] === undefined) data[category][k] = {};
+ if (data[category][k][v] === undefined) data[category][k][v] = 0;And increment the key for this key/value combination.
+ + data[category][k][item[k]]++;
+ }Increment the number of items classified
+ + total_count++;
+ };Generate a score of how well this item matches all +possible categories based on its attributes
+ + bayes_model.score = function(item) {Initialize an empty array of odds per category.
+ + var odds = {}, category;Iterate through each key in the item,
+then iterate through each category that has been used
+in previous calls to .train()
for (var k in item) {
+ var v = item[k];
+ for (category in data) {Create an empty object for storing key - value combinations +for this category.
+ + if (odds[category] === undefined) odds[category] = {};If this item doesn’t even have a property, it counts for nothing, +but if it does have the property that we’re looking for from +the item to categorize, it counts based on how popular it is +versus the whole population.
+ + if (data[category][k]) {
+ odds[category][k + '_' + v] = (data[category][k][v] || 0) / total_count;
+ } else {
+ odds[category][k + '_' + v] = 0;
+ }
+ }
+ }Set up a new object that will contain sums of these odds by category
+ + var odds_sums = {};
+
+ for (category in odds) {Tally all of the odds for each category-combination pair - +the non-existence of a category does not add anything to the +score.
+ + for (var combination in odds[category]) {
+ if (odds_sums[category] === undefined) odds_sums[category] = 0;
+ odds_sums[category] += odds[category][combination];
+ }
+ }
+
+ return odds_sums;
+ };Return the completed model.
+ + return bayes_model;
+ }is simply the result of adding all numbers +together, starting from zero.
+This runs on O(n), linear time in respect to the array
function sum(x) {
+ var value = 0;
+ for (var i = 0; i < x.length; i++) {
+ value += x[i];
+ }
+ return value;
+ }is the sum over the number of values
+This runs on O(n), linear time in respect to the array
function mean(x) {The mean of no numbers is null
+ + if (x.length === 0) return null;
+
+ return sum(x) / x.length;
+ }a mean function that is more useful for numbers in different +ranges.
+this is the nth root of the input numbers multiplied by each other
+This runs on O(n), linear time in respect to the array
function geometric_mean(x) {The mean of no numbers is null
+ + if (x.length === 0) return null;the starting value.
+ + var value = 1;
+
+ for (var i = 0; i < x.length; i++) {the geometric mean is only valid for positive numbers
+ + if (x[i] <= 0) return null;repeatedly multiply the value by each number
+ + value *= x[i];
+ }
+
+ return Math.pow(value, 1 / x.length);
+ }a mean function typically used to find the average of rates
+this is the reciprocal of the arithmetic mean of the reciprocals +of the input numbers
+This runs on O(n), linear time in respect to the array
function harmonic_mean(x) {The mean of no numbers is null
+ + if (x.length === 0) return null;
+
+ var reciprocal_sum = 0;
+
+ for (var i = 0; i < x.length; i++) {the harmonic mean is only valid for positive numbers
+ + if (x[i] <= 0) return null;
+
+ reciprocal_sum += 1 / x[i];
+ }divide n by the the reciprocal sum
+ + return x.length / reciprocal_sum;
+ }This is simply the minimum number in the set.
+This runs on O(n), linear time in respect to the array
function min(x) {
+ var value;
+ for (var i = 0; i < x.length; i++) {On the first iteration of this loop, min is +undefined and is thus made the minimum element in the array
+ + if (x[i] < value || value === undefined) value = x[i];
+ }
+ return value;
+ }This is simply the maximum number in the set.
+This runs on O(n), linear time in respect to the array
function max(x) {
+ var value;
+ for (var i = 0; i < x.length; i++) {On the first iteration of this loop, max is +undefined and is thus made the maximum element in the array
+ + if (x[i] > value || value === undefined) value = x[i];
+ }
+ return value;
+ } function variance(x) {The variance of no numbers is null
+ + if (x.length === 0) return null;
+
+ var mean_value = mean(x),
+ deviations = [];Make a list of squared deviations from the mean.
+ + for (var i = 0; i < x.length; i++) {
+ deviations.push(Math.pow(x[i] - mean_value, 2));
+ }Find the mean value of that list
+ + return mean(deviations);
+ } function standard_deviation(x) {The standard deviation of no numbers is null
+ + if (x.length === 0) return null;
+
+ return Math.sqrt(variance(x));
+ }The sum of deviations to the Nth power. +When n=2 it’s the sum of squared deviations. +When n=3 it’s the sum of cubed deviations.
+depends on mean()
function sum_nth_power_deviations(x, n) {
+ var mean_value = mean(x),
+ sum = 0;
+
+ for (var i = 0; i < x.length; i++) {
+ sum += Math.pow(x[i] - mean_value, n);
+ }
+
+ return sum;
+ }is the sum of squared deviations from the mean
+depends on sum_nth_power_deviations
function sample_variance(x) {The variance of no numbers is null
+ + if (x.length <= 1) return null;
+
+ var sum_squared_deviations_value = sum_nth_power_deviations(x, 2);Find the mean value of that list
+ + return sum_squared_deviations_value / (x.length - 1);
+ }is just the square root of the variance.
+depends on sample_variance()
function sample_standard_deviation(x) {The standard deviation of no numbers is null
+ + if (x.length <= 1) return null;
+
+ return Math.sqrt(sample_variance(x));
+ }sample covariance of two datasets: +how much do the two datasets move together? +x and y are two datasets, represented as arrays of numbers.
+depends on mean()
function sample_covariance(x, y) {The two datasets must have the same length which must be more than 1
+ + if (x.length <= 1 || x.length != y.length){
+ return null;
+ }determine the mean of each dataset so that we can judge each +value of the dataset fairly as the difference from the mean. this +way, if one dataset is [1, 2, 3] and [2, 3, 4], their covariance +does not suffer because of the difference in absolute values
+ + var xmean = mean(x),
+ ymean = mean(y),
+ sum = 0;for each pair of values, the covariance increases when their +difference from the mean is associated - if both are well above +or if both are well below +the mean, the covariance increases significantly.
+ + for (var i = 0; i < x.length; i++){
+ sum += (x[i] - xmean) * (y[i] - ymean);
+ }the covariance is weighted by the length of the datasets.
+ + return sum / (x.length - 1);
+ }Gets a measure of how correlated two datasets are, between -1 and 1
+depends on sample_standard_deviation() and sample_covariance()
function sample_correlation(x, y) {
+ var cov = sample_covariance(x, y),
+ xstd = sample_standard_deviation(x),
+ ystd = sample_standard_deviation(y);
+
+ if (cov === null || xstd === null || ystd === null) {
+ return null;
+ }
+
+ return cov / xstd / ystd;
+ }The middle number of a list. This is often a good indicator of ‘the middle’
+when there are outliers that skew the mean() value.
function median(x) {The median of an empty list is null
+ + if (x.length === 0) return null;Sorting the array makes it easy to find the center, but
+use .slice() to ensure the original array x is not modified
var sorted = x.slice().sort(function (a, b) { return a - b; });If the length of the list is odd, it’s the central number
+ + if (sorted.length % 2 === 1) {
+ return sorted[(sorted.length - 1) / 2];Otherwise, the median is the average of the two numbers +at the center of the list
+ + } else {
+ var a = sorted[(sorted.length / 2) - 1];
+ var b = sorted[(sorted.length / 2)];
+ return (a + b) / 2;
+ }
+ }The mode is the number that appears in a list the highest number of times. +There can be multiple modes in a list: in the event of a tie, this +algorithm will return the most recently seen mode.
+This implementation is inspired by science.js
+This runs on O(n), linear time in respect to the array
function mode(x) {Handle edge cases: +The median of an empty list is null
+ + if (x.length === 0) return null;
+ else if (x.length === 1) return x[0];Sorting the array lets us iterate through it below and be sure +that every time we see a new number it’s new and we’ll never +see the same number twice
+ + var sorted = x.slice().sort(function (a, b) { return a - b; });This assumes it is dealing with an array of size > 1, since size +0 and 1 are handled immediately. Hence it starts at index 1 in the +array.
+ + var last = sorted[0],store the mode as we find new modes
+ + value,store how many times we’ve seen the mode
+ + max_seen = 0,how many times the current candidate for the mode +has been seen
+ + seen_this = 1;end at sorted.length + 1 to fix the case in which the mode is +the highest number that occurs in the sequence. the last iteration +compares sorted[i], which is undefined, to the highest number +in the series
+ + for (var i = 1; i < sorted.length + 1; i++) {we’re seeing a new number pass by
+ + if (sorted[i] !== last) {the last number is the new mode since we saw it more +often than the old one
+ + if (seen_this > max_seen) {
+ max_seen = seen_this;
+ value = last;
+ }
+ seen_this = 1;
+ last = sorted[i];if this isn’t a new number, it’s one more occurrence of +the potential mode
+ + } else { seen_this++; }
+ }
+ return value;
+ }This is to compute a one-sample t-test, comparing the mean +of a sample to a known value, x.
+in this case, we’re trying to determine whether the
+population mean is equal to the value that we know, which is x
+here. usually the results here are used to look up a
+p-value, which, for
+a certain level of significance, will let you determine that the
+null hypothesis can or cannot be rejected.
Depends on standard_deviation() and mean()
function t_test(sample, x) {The mean of the sample
+ + var sample_mean = mean(sample);The standard deviation of the sample
+ + var sd = standard_deviation(sample);Square root the length of the sample
+ + var rootN = Math.sqrt(sample.length);Compute the known value against the sample, +returning the t value
+ + return (sample_mean - x) / (sd / rootN);
+ }This is to compute two sample t-test.
+Tests whether “mean(X)-mean(Y) = difference”, (
+in the most common case, we often have difference == 0 to test if two samples
+are likely to be taken from populations with the same mean value) with
+no prior knowledge on standard deviations of both samples
+other than the fact that they have the same standard deviation.
Usually the results here are used to look up a +p-value, which, for +a certain level of significance, will let you determine that the +null hypothesis can or cannot be rejected.
+diff can be omitted if it equals 0.
This is used to confirm or deny
+a null hypothesis that the two populations that have been sampled into
+sample_x and sample_y are equal to each other.
Depends on sample_variance() and mean()
function t_test_two_sample(sample_x, sample_y, difference) {
+ var n = sample_x.length,
+ m = sample_y.length;If either sample doesn’t actually have any values, we can’t
+compute this at all, so we return null.
if (!n || !m) return null ;default difference (mu) is zero
+ + if (!difference) difference = 0;
+
+ var meanX = mean(sample_x),
+ meanY = mean(sample_y);
+
+ var weightedVariance = ((n - 1) * sample_variance(sample_x) +
+ (m - 1) * sample_variance(sample_y)) / (n + m - 2);
+
+ return (meanX - meanY - difference) /
+ Math.sqrt(weightedVariance * (1 / n + 1 / m));
+ }Split an array into chunks of a specified size. This function +has the same behavior as PHP’s array_chunk +function, and thus will insert smaller-sized chunks at the end if +the input size is not divisible by the chunk size.
+sample is expected to be an array, and chunkSize a number.
+The sample array can contain any kind of data.
function chunk(sample, chunkSize) {a list of result chunks, as arrays in an array
+ + var output = [];chunkSize must be zero or higher - otherwise the loop below,
+in which we call start += chunkSize, will loop infinitely.
+So, we’ll detect and return null in that case to indicate
+invalid input.
if (chunkSize <= 0) {
+ return null;
+ }start is the index at which .slice will start selecting
+new array elements
for (var start = 0; start < sample.length; start += chunkSize) {for each chunk, slice that part of the array and add it
+to the output. The .slice function does not change
+the original array.
output.push(sample.slice(start, start + chunkSize));
+ }
+ return output;
+ }A Fisher-Yates shuffle +in-place - which means that it will change the order of the original +array by reference.
+ + function shuffle_in_place(sample, randomSource) {a custom random number source can be provided if you want to use +a fixed seed or another random number generator, like +random-js
+ + randomSource = randomSource || Math.random;store the current length of the sample to determine +when no elements remain to shuffle.
+ + var length = sample.length;temporary is used to hold an item when it is being +swapped between indices.
+ + var temporary;The index to swap at each stage.
+ + var index;While there are still items to shuffle
+ + while (length > 0) {chose a random index within the subset of the array +that is not yet shuffled
+ + index = Math.floor(randomSource() * length--);store the value that we’ll move temporarily
+ + temporary = sample[length];swap the value at sample[length] with sample[index]
sample[length] = sample[index];
+ sample[index] = temporary;
+ }
+
+ return sample;
+ }A Fisher-Yates shuffle +is a fast way to create a random permutation of a finite set.
+ + function shuffle(sample, randomSource) {slice the original array so that it is not modified
+ + sample = sample.slice();and then shuffle that shallow-copied array, in place
+ + return shuffle_in_place(sample.slice(), randomSource);
+ } function sample(array, n, randomSource) {shuffle the original array using a fisher-yates shuffle
+ + var shuffled = shuffle(array, randomSource);and then return a subset of it - the first n elements.
return shuffled.slice(0, n);
+ }This is a population quantile, since we assume to know the entire +dataset in this library. Thus I’m trying to follow the +Quantiles of a Population +algorithm from wikipedia.
+Sample is a one-dimensional array of numbers, +and p is either a decimal number from 0 to 1 or an array of decimal +numbers from 0 to 1. +In terms of a k/q quantile, p = k/q - it’s just dealing with fractions or dealing +with decimal values. +When p is an array, the result of the function is also an array containing the appropriate +quantiles in input order
+ + function quantile(sample, p) {We can’t derive quantiles from an empty list
+ + if (sample.length === 0) return null;Sort a copy of the array. We’ll need a sorted array to index +the values in sorted order.
+ + var sorted = sample.slice().sort(function (a, b) { return a - b; });
+
+ if (p.length) {Initialize the result array
+ + var results = [];For each requested quantile
+ + for (var i = 0; i < p.length; i++) {
+ results[i] = quantile_sorted(sorted, p[i]);
+ }
+ return results;
+ } else {
+ return quantile_sorted(sorted, p);
+ }
+ }This is the internal implementation of quantiles: when you know +that the order is sorted, you don’t need to re-sort it, and the computations +are much faster.
+ + function quantile_sorted(sample, p) {
+ var idx = (sample.length) * p;
+ if (p < 0 || p > 1) {
+ return null;
+ } else if (p === 1) {If p is 1, directly return the last element
+ + return sample[sample.length - 1];
+ } else if (p === 0) {If p is 0, directly return the first element
+ + return sample[0];
+ } else if (idx % 1 !== 0) {If p is not integer, return the next element in array
+ + return sample[Math.ceil(idx) - 1];
+ } else if (sample.length % 2 === 0) {If the list has even-length, we’ll take the average of this number +and the next value, if there is one
+ + return (sample[idx - 1] + sample[idx]) / 2;
+ } else {Finally, in the simple case of an integer value +with an odd-length list, return the sample value at the index.
+ + return sample[idx];
+ }
+ }A measure of statistical dispersion, or how scattered, spread, or +concentrated a distribution is. It’s computed as the difference between +the third quartile and first quartile.
+ + function iqr(sample) {We can’t derive quantiles from an empty list
+ + if (sample.length === 0) return null;Interquartile range is the span between the upper quartile,
+at 0.75, and lower quartile, 0.25
return quantile(sample, 0.75) - quantile(sample, 0.25);
+ }The Median Absolute Deviation (MAD) is a robust measure of statistical +dispersion. It is more resilient to outliers than the standard deviation.
+ + function mad(x) {The mad of nothing is null
+ + if (!x || x.length === 0) return null;
+
+ var median_value = median(x),
+ median_absolute_deviations = [];Make a list of absolute deviations from the median
+ + for (var i = 0; i < x.length; i++) {
+ median_absolute_deviations.push(Math.abs(x[i] - median_value));
+ }Find the median value of that list
+ + return median(median_absolute_deviations);
+ }Compute the matrices required for Jenks breaks. These matrices
+can be used for any classing of data with classes <= n_classes
function jenksMatrices(data, n_classes) {in the original implementation, these matrices are referred to
+as LC and OP
var lower_class_limits = [],
+ variance_combinations = [],loop counters
+ + i, j,the variance, as computed at each step in the calculation
+ + variance = 0;Initialize and fill each matrix with zeroes
+ + for (i = 0; i < data.length + 1; i++) {
+ var tmp1 = [], tmp2 = [];despite these arrays having the same values, we need +to keep them separate so that changing one does not change +the other
+ + for (j = 0; j < n_classes + 1; j++) {
+ tmp1.push(0);
+ tmp2.push(0);
+ }
+ lower_class_limits.push(tmp1);
+ variance_combinations.push(tmp2);
+ }
+
+ for (i = 1; i < n_classes + 1; i++) {
+ lower_class_limits[1][i] = 1;
+ variance_combinations[1][i] = 0;in the original implementation, 9999999 is used but
+since Javascript has Infinity, we use that.
for (j = 2; j < data.length + 1; j++) {
+ variance_combinations[j][i] = Infinity;
+ }
+ }
+
+ for (var l = 2; l < data.length + 1; l++) {SZ originally. this is the sum of the values seen thus
+far when calculating variance.
var sum = 0,ZSQ originally. the sum of squares of values seen
+thus far
sum_squares = 0,WT originally. This is the number of
w = 0,IV originally
i4 = 0;in several instances, you could say Math.pow(x, 2)
+instead of x * x, but this is slower in some browsers
+introduces an unnecessary concept.
for (var m = 1; m < l + 1; m++) {III originally
var lower_class_limit = l - m + 1,
+ val = data[lower_class_limit - 1];here we’re estimating variance for each potential classing
+of the data, for each potential number of classes. w
+is the number of data points considered so far.
w++;increase the current sum and sum-of-squares
+ + sum += val;
+ sum_squares += val * val;the variance at this point in the sequence is the difference +between the sum of squares and the total x 2, over the number +of samples.
+ + variance = sum_squares - (sum * sum) / w;
+
+ i4 = lower_class_limit - 1;
+
+ if (i4 !== 0) {
+ for (j = 2; j < n_classes + 1; j++) {if adding this element to an existing class
+will increase its variance beyond the limit, break
+the class at this point, setting the lower_class_limit
+at this point.
if (variance_combinations[l][j] >=
+ (variance + variance_combinations[i4][j - 1])) {
+ lower_class_limits[l][j] = lower_class_limit;
+ variance_combinations[l][j] = variance +
+ variance_combinations[i4][j - 1];
+ }
+ }
+ }
+ }
+
+ lower_class_limits[l][1] = 1;
+ variance_combinations[l][1] = variance;
+ }return the two matrices. for just providing breaks, only
+lower_class_limits is needed, but variances can be useful to
+evaluate goodness of fit.
return {
+ lower_class_limits: lower_class_limits,
+ variance_combinations: variance_combinations
+ };
+ }the second part of the jenks recipe: take the calculated matrices +and derive an array of n breaks.
+ + function jenksBreaks(data, lower_class_limits, n_classes) {
+
+ var k = data.length - 1,
+ kclass = [],
+ countNum = n_classes;the calculation of classes will never include the upper and +lower bounds, so we need to explicitly set them
+ + kclass[n_classes] = data[data.length - 1];
+ kclass[0] = data[0];the lower_class_limits matrix is used as indices into itself
+here: the k variable is reused in each iteration.
while (countNum > 1) {
+ kclass[countNum - 1] = data[lower_class_limits[k][countNum] - 2];
+ k = lower_class_limits[k][countNum] - 1;
+ countNum--;
+ }
+
+ return kclass;
+ }Implementations: 1 (python), +2 (buggy), +3 (works)
+Depends on jenksBreaks() and jenksMatrices()
function jenks(data, n_classes) {
+
+ if (n_classes > data.length) return null;sort data in numerical order, since this is expected +by the matrices function
+ + data = data.slice().sort(function (a, b) { return a - b; });get our basic matrices
+ + var matrices = jenksMatrices(data, n_classes),we only need lower class limits here
+ + lower_class_limits = matrices.lower_class_limits;extract n_classes out of the computed matrices
+ + return jenksBreaks(data, lower_class_limits, n_classes);
+
+ }A measure of the extent to which a probability distribution of a +real-valued random variable “leans” to one side of the mean. +The skewness value can be positive or negative, or even undefined.
+Implementation is based on the adjusted Fisher-Pearson standardized +moment coefficient, which is the version found in Excel and several +statistical packages including Minitab, SAS and SPSS.
+Depends on sum_nth_power_deviations() and sample_standard_deviation
function sample_skewness(x) {The skewness of less than three arguments is null
+ + if (x.length < 3) return null;
+
+ var n = x.length,
+ cubed_s = Math.pow(sample_standard_deviation(x), 3),
+ sum_cubed_deviations = sum_nth_power_deviations(x, 3);
+
+ return n * sum_cubed_deviations / ((n - 1) * (n - 2) * cubed_s);
+ }A standard normal table, also called the unit normal table or Z table, +is a mathematical table for the values of Φ (phi), which are the values of +the cumulative distribution function of the normal distribution. +It is used to find the probability that a statistic is observed below, +above, or between values on the standard normal distribution, and by +extension, any normal distribution.
+The probabilities are taken from http://en.wikipedia.org/wiki/Standard_normal_table +The table used is the cumulative, and not cumulative from 0 to mean +(even though the latter has 5 digits precision, instead of 4).
+ + var standard_normal_table = [
+ /* z 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 */
+ /* 0.0 */
+ 0.5000, 0.5040, 0.5080, 0.5120, 0.5160, 0.5199, 0.5239, 0.5279, 0.5319, 0.5359,
+ /* 0.1 */
+ 0.5398, 0.5438, 0.5478, 0.5517, 0.5557, 0.5596, 0.5636, 0.5675, 0.5714, 0.5753,
+ /* 0.2 */
+ 0.5793, 0.5832, 0.5871, 0.5910, 0.5948, 0.5987, 0.6026, 0.6064, 0.6103, 0.6141,
+ /* 0.3 */
+ 0.6179, 0.6217, 0.6255, 0.6293, 0.6331, 0.6368, 0.6406, 0.6443, 0.6480, 0.6517,
+ /* 0.4 */
+ 0.6554, 0.6591, 0.6628, 0.6664, 0.6700, 0.6736, 0.6772, 0.6808, 0.6844, 0.6879,
+ /* 0.5 */
+ 0.6915, 0.6950, 0.6985, 0.7019, 0.7054, 0.7088, 0.7123, 0.7157, 0.7190, 0.7224,
+ /* 0.6 */
+ 0.7257, 0.7291, 0.7324, 0.7357, 0.7389, 0.7422, 0.7454, 0.7486, 0.7517, 0.7549,
+ /* 0.7 */
+ 0.7580, 0.7611, 0.7642, 0.7673, 0.7704, 0.7734, 0.7764, 0.7794, 0.7823, 0.7852,
+ /* 0.8 */
+ 0.7881, 0.7910, 0.7939, 0.7967, 0.7995, 0.8023, 0.8051, 0.8078, 0.8106, 0.8133,
+ /* 0.9 */
+ 0.8159, 0.8186, 0.8212, 0.8238, 0.8264, 0.8289, 0.8315, 0.8340, 0.8365, 0.8389,
+ /* 1.0 */
+ 0.8413, 0.8438, 0.8461, 0.8485, 0.8508, 0.8531, 0.8554, 0.8577, 0.8599, 0.8621,
+ /* 1.1 */
+ 0.8643, 0.8665, 0.8686, 0.8708, 0.8729, 0.8749, 0.8770, 0.8790, 0.8810, 0.8830,
+ /* 1.2 */
+ 0.8849, 0.8869, 0.8888, 0.8907, 0.8925, 0.8944, 0.8962, 0.8980, 0.8997, 0.9015,
+ /* 1.3 */
+ 0.9032, 0.9049, 0.9066, 0.9082, 0.9099, 0.9115, 0.9131, 0.9147, 0.9162, 0.9177,
+ /* 1.4 */
+ 0.9192, 0.9207, 0.9222, 0.9236, 0.9251, 0.9265, 0.9279, 0.9292, 0.9306, 0.9319,
+ /* 1.5 */
+ 0.9332, 0.9345, 0.9357, 0.9370, 0.9382, 0.9394, 0.9406, 0.9418, 0.9429, 0.9441,
+ /* 1.6 */
+ 0.9452, 0.9463, 0.9474, 0.9484, 0.9495, 0.9505, 0.9515, 0.9525, 0.9535, 0.9545,
+ /* 1.7 */
+ 0.9554, 0.9564, 0.9573, 0.9582, 0.9591, 0.9599, 0.9608, 0.9616, 0.9625, 0.9633,
+ /* 1.8 */
+ 0.9641, 0.9649, 0.9656, 0.9664, 0.9671, 0.9678, 0.9686, 0.9693, 0.9699, 0.9706,
+ /* 1.9 */
+ 0.9713, 0.9719, 0.9726, 0.9732, 0.9738, 0.9744, 0.9750, 0.9756, 0.9761, 0.9767,
+ /* 2.0 */
+ 0.9772, 0.9778, 0.9783, 0.9788, 0.9793, 0.9798, 0.9803, 0.9808, 0.9812, 0.9817,
+ /* 2.1 */
+ 0.9821, 0.9826, 0.9830, 0.9834, 0.9838, 0.9842, 0.9846, 0.9850, 0.9854, 0.9857,
+ /* 2.2 */
+ 0.9861, 0.9864, 0.9868, 0.9871, 0.9875, 0.9878, 0.9881, 0.9884, 0.9887, 0.9890,
+ /* 2.3 */
+ 0.9893, 0.9896, 0.9898, 0.9901, 0.9904, 0.9906, 0.9909, 0.9911, 0.9913, 0.9916,
+ /* 2.4 */
+ 0.9918, 0.9920, 0.9922, 0.9925, 0.9927, 0.9929, 0.9931, 0.9932, 0.9934, 0.9936,
+ /* 2.5 */
+ 0.9938, 0.9940, 0.9941, 0.9943, 0.9945, 0.9946, 0.9948, 0.9949, 0.9951, 0.9952,
+ /* 2.6 */
+ 0.9953, 0.9955, 0.9956, 0.9957, 0.9959, 0.9960, 0.9961, 0.9962, 0.9963, 0.9964,
+ /* 2.7 */
+ 0.9965, 0.9966, 0.9967, 0.9968, 0.9969, 0.9970, 0.9971, 0.9972, 0.9973, 0.9974,
+ /* 2.8 */
+ 0.9974, 0.9975, 0.9976, 0.9977, 0.9977, 0.9978, 0.9979, 0.9979, 0.9980, 0.9981,
+ /* 2.9 */
+ 0.9981, 0.9982, 0.9982, 0.9983, 0.9984, 0.9984, 0.9985, 0.9985, 0.9986, 0.9986,
+ /* 3.0 */
+ 0.9987, 0.9987, 0.9987, 0.9988, 0.9988, 0.9989, 0.9989, 0.9989, 0.9990, 0.9990
+ ];Since probability tables cannot be +printed for every normal distribution, as there are an infinite variety +of normal distributions, it is common practice to convert a normal to a +standard normal and then use the standard normal table to find probabilities
+ + function cumulative_std_normal_probability(z) {Calculate the position of this value.
+ + var absZ = Math.abs(z),Each row begins with a different +significant digit: 0.5, 0.6, 0.7, and so on. So the row is simply +this value’s significant digit: 0.567 will be in row 0, so row=0, +0.643 will be in row 1, so row=10.
+ + row = Math.floor(absZ * 10),
+ column = 10 * (Math.floor(absZ * 100) / 10 - Math.floor(absZ * 100 / 10)),
+ index = Math.min((row * 10) + column, standard_normal_table.length - 1);The index we calculate must be in the table as a positive value, +but we still pay attention to whether the input is positive +or negative, and flip the output value as a last step.
+ + if (z >= 0) {
+ return standard_normal_table[index];
+ } else {due to floating-point arithmetic, values in the table with +4 significant figures can nevertheless end up as repeating +fractions when they’re computed here.
+ + return +(1 - standard_normal_table[index]).toFixed(4);
+ }
+ }The standard score is the number of standard deviations an observation +or datum is above or below the mean. Thus, a positive standard score +represents a datum above the mean, while a negative standard score +represents a datum below the mean. It is a dimensionless quantity +obtained by subtracting the population mean from an individual raw +score and then dividing the difference by the population standard +deviation.
+The z-score is only defined if one knows the population parameters; +if one only has a sample set, then the analogous computation with +sample mean and sample standard deviation yields the +Student’s t-statistic.
+ + function z_score(x, mean, standard_deviation) {
+ return (x - mean) / standard_deviation;
+ }We use ε, epsilon, as a stopping criterion when we want to iterate
+until we’re “close enough”.
var epsilon = 0.0001;A factorial, usually written n!, is the product of all positive +integers less than or equal to n. Often factorial is implemented +recursively, but this iterative approach is significantly faster +and simpler.
+ + function factorial(n) {factorial is mathematically undefined for negative numbers
+ + if (n < 0 ) { return null; }typically you’ll expand the factorial function going down, like +5! = 5 4 3 2 1. This is going in the opposite direction, +counting from 2 up to the number in question, and since anything +multiplied by 1 is itself, the loop only needs to start at 2.
+ + var accumulator = 1;
+ for (var i = 2; i <= n; i++) {for each number up to and including the number n, multiply
+the accumulator my that number.
accumulator *= i;
+ }
+ return accumulator;
+ }The Bernoulli distribution
+is the probability discrete
+distribution of a random variable which takes value 1 with success
+probability p and value 0 with failure
+probability q = 1 - p. It can be used, for example, to represent the
+toss of a coin, where “1” is defined to mean “heads” and “0” is defined
+to mean “tails” (or vice versa). It is
+a special case of a Binomial Distribution
+where n = 1.
function bernoulli_distribution(p) {Check that p is a valid probability (0 ≤ p ≤ 1)
if (p < 0 || p > 1 ) { return null; }
+
+ return binomial_distribution(1, p);
+ }The Binomial Distribution is the discrete probability
+distribution of the number of successes in a sequence of n independent yes/no experiments, each of which yields
+success with probability probability. Such a success/failure experiment is also called a Bernoulli experiment or
+Bernoulli trial; when trials = 1, the Binomial Distribution is a Bernoulli Distribution.
function binomial_distribution(trials, probability) {Check that p is a valid probability (0 ≤ p ≤ 1),
+that n is an integer, strictly positive.
if (probability < 0 || probability > 1 ||
+ trials <= 0 || trials % 1 !== 0) {
+ return null;
+ } function probability_mass(x, trials, probability) {
+ return factorial(trials) /
+ (factorial(x) * factorial(trials - x)) *
+ (Math.pow(probability, x) * Math.pow(1 - probability, trials - x));
+ }We initialize x, the random variable, and accumulator, an accumulator
+for the cumulative distribution function to 0. distribution_functions
+is the object we’ll return with the probability_of_x and the
+cumulative_probability_of_x, as well as the calculated mean &
+variance. We iterate until the cumulative_probability_of_x is
+within epsilon of 1.0.
var x = 0,
+ cumulative_probability = 0,
+ cells = {};This algorithm iterates through each potential outcome,
+until the cumulative_probability is very close to 1, at
+which point we’ve defined the vast majority of outcomes
do {
+ cells[x] = probability_mass(x, trials, probability);
+ cumulative_probability += cells[x];
+ x++;when the cumulative_probability is nearly 1, we’ve calculated +the useful range of this distribution
+ + } while (cumulative_probability < 1 - epsilon);
+
+ return cells;
+ }The Poisson Distribution +is a discrete probability distribution that expresses the probability +of a given number of events occurring in a fixed interval of time +and/or space if these events occur with a known average rate and +independently of the time since the last event.
+The Poisson Distribution is characterized by the strictly positive
+mean arrival or occurrence rate, λ.
function poisson_distribution(lambda) {Check that lambda is strictly positive
+ + if (lambda <= 0) { return null; }our current place in the distribution
+ + var x = 0,and we keep track of the current cumulative probability, in +order to know when to stop calculating chances.
+ + cumulative_probability = 0,the calculated cells to be returned
+ + cells = {}; function probability_mass(x, lambda) {
+ return (Math.pow(Math.E, -lambda) * Math.pow(lambda, x)) /
+ factorial(x);
+ }This algorithm iterates through each potential outcome,
+until the cumulative_probability is very close to 1, at
+which point we’ve defined the vast majority of outcomes
do {
+ cells[x] = probability_mass(x, lambda);
+ cumulative_probability += cells[x];
+ x++;when the cumulative_probability is nearly 1, we’ve calculated +the useful range of this distribution
+ + } while (cumulative_probability < 1 - epsilon);
+
+ return cells;
+ }The χ2 (Chi-Squared) Distribution is used in the common +chi-squared tests for goodness of fit of an observed distribution to a theoretical one, the independence of two +criteria of classification of qualitative data, and in confidence interval estimation for a population standard +deviation of a normal distribution from a sample standard deviation.
+Values from Appendix 1, Table III of William W. Hines & Douglas C. Montgomery, “Probability and Statistics in +Engineering and Management Science”, Wiley (1980).
+ + var chi_squared_distribution_table = {
+ 1: { 0.995: 0.00, 0.99: 0.00, 0.975: 0.00, 0.95: 0.00, 0.9: 0.02, 0.5: 0.45, 0.1: 2.71, 0.05: 3.84, 0.025: 5.02, 0.01: 6.63, 0.005: 7.88 },
+ 2: { 0.995: 0.01, 0.99: 0.02, 0.975: 0.05, 0.95: 0.10, 0.9: 0.21, 0.5: 1.39, 0.1: 4.61, 0.05: 5.99, 0.025: 7.38, 0.01: 9.21, 0.005: 10.60 },
+ 3: { 0.995: 0.07, 0.99: 0.11, 0.975: 0.22, 0.95: 0.35, 0.9: 0.58, 0.5: 2.37, 0.1: 6.25, 0.05: 7.81, 0.025: 9.35, 0.01: 11.34, 0.005: 12.84 },
+ 4: { 0.995: 0.21, 0.99: 0.30, 0.975: 0.48, 0.95: 0.71, 0.9: 1.06, 0.5: 3.36, 0.1: 7.78, 0.05: 9.49, 0.025: 11.14, 0.01: 13.28, 0.005: 14.86 },
+ 5: { 0.995: 0.41, 0.99: 0.55, 0.975: 0.83, 0.95: 1.15, 0.9: 1.61, 0.5: 4.35, 0.1: 9.24, 0.05: 11.07, 0.025: 12.83, 0.01: 15.09, 0.005: 16.75 },
+ 6: { 0.995: 0.68, 0.99: 0.87, 0.975: 1.24, 0.95: 1.64, 0.9: 2.20, 0.5: 5.35, 0.1: 10.65, 0.05: 12.59, 0.025: 14.45, 0.01: 16.81, 0.005: 18.55 },
+ 7: { 0.995: 0.99, 0.99: 1.25, 0.975: 1.69, 0.95: 2.17, 0.9: 2.83, 0.5: 6.35, 0.1: 12.02, 0.05: 14.07, 0.025: 16.01, 0.01: 18.48, 0.005: 20.28 },
+ 8: { 0.995: 1.34, 0.99: 1.65, 0.975: 2.18, 0.95: 2.73, 0.9: 3.49, 0.5: 7.34, 0.1: 13.36, 0.05: 15.51, 0.025: 17.53, 0.01: 20.09, 0.005: 21.96 },
+ 9: { 0.995: 1.73, 0.99: 2.09, 0.975: 2.70, 0.95: 3.33, 0.9: 4.17, 0.5: 8.34, 0.1: 14.68, 0.05: 16.92, 0.025: 19.02, 0.01: 21.67, 0.005: 23.59 },
+ 10: { 0.995: 2.16, 0.99: 2.56, 0.975: 3.25, 0.95: 3.94, 0.9: 4.87, 0.5: 9.34, 0.1: 15.99, 0.05: 18.31, 0.025: 20.48, 0.01: 23.21, 0.005: 25.19 },
+ 11: { 0.995: 2.60, 0.99: 3.05, 0.975: 3.82, 0.95: 4.57, 0.9: 5.58, 0.5: 10.34, 0.1: 17.28, 0.05: 19.68, 0.025: 21.92, 0.01: 24.72, 0.005: 26.76 },
+ 12: { 0.995: 3.07, 0.99: 3.57, 0.975: 4.40, 0.95: 5.23, 0.9: 6.30, 0.5: 11.34, 0.1: 18.55, 0.05: 21.03, 0.025: 23.34, 0.01: 26.22, 0.005: 28.30 },
+ 13: { 0.995: 3.57, 0.99: 4.11, 0.975: 5.01, 0.95: 5.89, 0.9: 7.04, 0.5: 12.34, 0.1: 19.81, 0.05: 22.36, 0.025: 24.74, 0.01: 27.69, 0.005: 29.82 },
+ 14: { 0.995: 4.07, 0.99: 4.66, 0.975: 5.63, 0.95: 6.57, 0.9: 7.79, 0.5: 13.34, 0.1: 21.06, 0.05: 23.68, 0.025: 26.12, 0.01: 29.14, 0.005: 31.32 },
+ 15: { 0.995: 4.60, 0.99: 5.23, 0.975: 6.27, 0.95: 7.26, 0.9: 8.55, 0.5: 14.34, 0.1: 22.31, 0.05: 25.00, 0.025: 27.49, 0.01: 30.58, 0.005: 32.80 },
+ 16: { 0.995: 5.14, 0.99: 5.81, 0.975: 6.91, 0.95: 7.96, 0.9: 9.31, 0.5: 15.34, 0.1: 23.54, 0.05: 26.30, 0.025: 28.85, 0.01: 32.00, 0.005: 34.27 },
+ 17: { 0.995: 5.70, 0.99: 6.41, 0.975: 7.56, 0.95: 8.67, 0.9: 10.09, 0.5: 16.34, 0.1: 24.77, 0.05: 27.59, 0.025: 30.19, 0.01: 33.41, 0.005: 35.72 },
+ 18: { 0.995: 6.26, 0.99: 7.01, 0.975: 8.23, 0.95: 9.39, 0.9: 10.87, 0.5: 17.34, 0.1: 25.99, 0.05: 28.87, 0.025: 31.53, 0.01: 34.81, 0.005: 37.16 },
+ 19: { 0.995: 6.84, 0.99: 7.63, 0.975: 8.91, 0.95: 10.12, 0.9: 11.65, 0.5: 18.34, 0.1: 27.20, 0.05: 30.14, 0.025: 32.85, 0.01: 36.19, 0.005: 38.58 },
+ 20: { 0.995: 7.43, 0.99: 8.26, 0.975: 9.59, 0.95: 10.85, 0.9: 12.44, 0.5: 19.34, 0.1: 28.41, 0.05: 31.41, 0.025: 34.17, 0.01: 37.57, 0.005: 40.00 },
+ 21: { 0.995: 8.03, 0.99: 8.90, 0.975: 10.28, 0.95: 11.59, 0.9: 13.24, 0.5: 20.34, 0.1: 29.62, 0.05: 32.67, 0.025: 35.48, 0.01: 38.93, 0.005: 41.40 },
+ 22: { 0.995: 8.64, 0.99: 9.54, 0.975: 10.98, 0.95: 12.34, 0.9: 14.04, 0.5: 21.34, 0.1: 30.81, 0.05: 33.92, 0.025: 36.78, 0.01: 40.29, 0.005: 42.80 },
+ 23: { 0.995: 9.26, 0.99: 10.20, 0.975: 11.69, 0.95: 13.09, 0.9: 14.85, 0.5: 22.34, 0.1: 32.01, 0.05: 35.17, 0.025: 38.08, 0.01: 41.64, 0.005: 44.18 },
+ 24: { 0.995: 9.89, 0.99: 10.86, 0.975: 12.40, 0.95: 13.85, 0.9: 15.66, 0.5: 23.34, 0.1: 33.20, 0.05: 36.42, 0.025: 39.36, 0.01: 42.98, 0.005: 45.56 },
+ 25: { 0.995: 10.52, 0.99: 11.52, 0.975: 13.12, 0.95: 14.61, 0.9: 16.47, 0.5: 24.34, 0.1: 34.28, 0.05: 37.65, 0.025: 40.65, 0.01: 44.31, 0.005: 46.93 },
+ 26: { 0.995: 11.16, 0.99: 12.20, 0.975: 13.84, 0.95: 15.38, 0.9: 17.29, 0.5: 25.34, 0.1: 35.56, 0.05: 38.89, 0.025: 41.92, 0.01: 45.64, 0.005: 48.29 },
+ 27: { 0.995: 11.81, 0.99: 12.88, 0.975: 14.57, 0.95: 16.15, 0.9: 18.11, 0.5: 26.34, 0.1: 36.74, 0.05: 40.11, 0.025: 43.19, 0.01: 46.96, 0.005: 49.65 },
+ 28: { 0.995: 12.46, 0.99: 13.57, 0.975: 15.31, 0.95: 16.93, 0.9: 18.94, 0.5: 27.34, 0.1: 37.92, 0.05: 41.34, 0.025: 44.46, 0.01: 48.28, 0.005: 50.99 },
+ 29: { 0.995: 13.12, 0.99: 14.26, 0.975: 16.05, 0.95: 17.71, 0.9: 19.77, 0.5: 28.34, 0.1: 39.09, 0.05: 42.56, 0.025: 45.72, 0.01: 49.59, 0.005: 52.34 },
+ 30: { 0.995: 13.79, 0.99: 14.95, 0.975: 16.79, 0.95: 18.49, 0.9: 20.60, 0.5: 29.34, 0.1: 40.26, 0.05: 43.77, 0.025: 46.98, 0.01: 50.89, 0.005: 53.67 },
+ 40: { 0.995: 20.71, 0.99: 22.16, 0.975: 24.43, 0.95: 26.51, 0.9: 29.05, 0.5: 39.34, 0.1: 51.81, 0.05: 55.76, 0.025: 59.34, 0.01: 63.69, 0.005: 66.77 },
+ 50: { 0.995: 27.99, 0.99: 29.71, 0.975: 32.36, 0.95: 34.76, 0.9: 37.69, 0.5: 49.33, 0.1: 63.17, 0.05: 67.50, 0.025: 71.42, 0.01: 76.15, 0.005: 79.49 },
+ 60: { 0.995: 35.53, 0.99: 37.48, 0.975: 40.48, 0.95: 43.19, 0.9: 46.46, 0.5: 59.33, 0.1: 74.40, 0.05: 79.08, 0.025: 83.30, 0.01: 88.38, 0.005: 91.95 },
+ 70: { 0.995: 43.28, 0.99: 45.44, 0.975: 48.76, 0.95: 51.74, 0.9: 55.33, 0.5: 69.33, 0.1: 85.53, 0.05: 90.53, 0.025: 95.02, 0.01: 100.42, 0.005: 104.22 },
+ 80: { 0.995: 51.17, 0.99: 53.54, 0.975: 57.15, 0.95: 60.39, 0.9: 64.28, 0.5: 79.33, 0.1: 96.58, 0.05: 101.88, 0.025: 106.63, 0.01: 112.33, 0.005: 116.32 },
+ 90: { 0.995: 59.20, 0.99: 61.75, 0.975: 65.65, 0.95: 69.13, 0.9: 73.29, 0.5: 89.33, 0.1: 107.57, 0.05: 113.14, 0.025: 118.14, 0.01: 124.12, 0.005: 128.30 },
+ 100: { 0.995: 67.33, 0.99: 70.06, 0.975: 74.22, 0.95: 77.93, 0.9: 82.36, 0.5: 99.33, 0.1: 118.50, 0.05: 124.34, 0.025: 129.56, 0.01: 135.81, 0.005: 140.17 }
+ };The χ2 (Chi-Squared) Goodness-of-Fit Test
+uses a measure of goodness of fit which is the sum of differences between observed and expected outcome frequencies
+(that is, counts of observations), each squared and divided by the number of observations expected given the
+hypothesized distribution. The resulting χ2 statistic, chi_squared, can be compared to the chi-squared distribution
+to determine the goodness of fit. In order to determine the degrees of freedom of the chi-squared distribution, one
+takes the total number of observed frequencies and subtracts the number of estimated parameters. The test statistic
+follows, approximately, a chi-square distribution with (k − c) degrees of freedom where k is the number of non-empty
+cells and c is the number of estimated parameters for the distribution.
function chi_squared_goodness_of_fit(data, distribution_type, significance) {Estimate from the sample data, a weighted mean.
+ + var input_mean = mean(data),Calculated value of the χ2 statistic.
+ + chi_squared = 0,Degrees of freedom, calculated as (number of class intervals - +number of hypothesized distribution parameters estimated - 1)
+ + degrees_of_freedom,Number of hypothesized distribution parameters estimated, expected to be supplied in the distribution test.
+Lose one degree of freedom for estimating lambda from the sample data.
c = 1,The hypothesized distribution. +Generate the hypothesized distribution.
+ + hypothesized_distribution = distribution_type(input_mean),
+ observed_frequencies = [],
+ expected_frequencies = [],
+ k;Create an array holding a histogram from the sample data, of
+the form { value: numberOfOcurrences }
for (var i = 0; i < data.length; i++) {
+ if (observed_frequencies[data[i]] === undefined) {
+ observed_frequencies[data[i]] = 0;
+ }
+ observed_frequencies[data[i]]++;
+ }The histogram we created might be sparse - there might be gaps +between values. So we iterate through the histogram, making +sure that instead of undefined, gaps have 0 values.
+ + for (i = 0; i < observed_frequencies.length; i++) {
+ if (observed_frequencies[i] === undefined) {
+ observed_frequencies[i] = 0;
+ }
+ }Create an array holding a histogram of expected data given the +sample size and hypothesized distribution.
+ + for (k in hypothesized_distribution) {
+ if (k in observed_frequencies) {
+ expected_frequencies[k] = hypothesized_distribution[k] * data.length;
+ }
+ }Working backward through the expected frequencies, collapse classes +if less than three observations are expected for a class. +This transformation is applied to the observed frequencies as well.
+ + for (k = expected_frequencies.length - 1; k >= 0; k--) {
+ if (expected_frequencies[k] < 3) {
+ expected_frequencies[k - 1] += expected_frequencies[k];
+ expected_frequencies.pop();
+
+ observed_frequencies[k - 1] += observed_frequencies[k];
+ observed_frequencies.pop();
+ }
+ }Iterate through the squared differences between observed & expected
+frequencies, accumulating the chi_squared statistic.
for (k = 0; k < observed_frequencies.length; k++) {
+ chi_squared += Math.pow(
+ observed_frequencies[k] - expected_frequencies[k], 2) /
+ expected_frequencies[k];
+ }Calculate degrees of freedom for this test and look it up in the
+chi_squared_distribution_table in order to
+accept or reject the goodness-of-fit of the hypothesized distribution.
degrees_of_freedom = observed_frequencies.length - c - 1;
+ return chi_squared_distribution_table[degrees_of_freedom][significance] < chi_squared;
+ }Mixin simple_statistics to a single Array instance if provided +or the Array native object if not. This is an optional +feature that lets you treat simple_statistics as a native feature +of Javascript.
+ + function mixin(array) {
+ var support = !!(Object.defineProperty && Object.defineProperties);
+ if (!support) throw new Error('without defineProperty, simple-statistics cannot be mixed in');only methods which work on basic arrays in a single step +are supported
+ + var arrayMethods = ['median', 'standard_deviation', 'sum',
+ 'sample_skewness',
+ 'mean', 'min', 'max', 'quantile', 'geometric_mean',
+ 'harmonic_mean'];create a closure with a method name so that a reference
+like arrayMethods[i] doesn’t follow the loop increment
function wrap(method) {
+ return function() {cast any arguments into an array, since they’re +natively objects
+ + var args = Array.prototype.slice.apply(arguments);make the first argument the array itself
+ + args.unshift(this);return the result of the ss method
+ + return ss[method].apply(ss, args);
+ };
+ }select object to extend
+ + var extending;
+ if (array) {create a shallow copy of the array so that our internal +operations do not change it by reference
+ + extending = array.slice();
+ } else {
+ extending = Array.prototype;
+ }for each array function, define a function that gets
+the array as the first argument.
+We use defineProperty
+because it allows these properties to be non-enumerable:
+for (var in x) loops will not run into problems with this
+implementation.
for (var i = 0; i < arrayMethods.length; i++) {
+ Object.defineProperty(extending, arrayMethods[i], {
+ value: wrap(arrayMethods[i]),
+ configurable: true,
+ enumerable: false,
+ writable: true
+ });
+ }
+
+ return extending;
+ }
+
+ ss.linear_regression = linear_regression;
+ ss.standard_deviation = standard_deviation;
+ ss.r_squared = r_squared;
+ ss.median = median;
+ ss.mean = mean;
+ ss.mode = mode;
+ ss.min = min;
+ ss.max = max;
+ ss.sum = sum;
+ ss.quantile = quantile;
+ ss.quantile_sorted = quantile_sorted;
+ ss.iqr = iqr;
+ ss.mad = mad;
+
+ ss.chunk = chunk;
+ ss.shuffle = shuffle;
+ ss.shuffle_in_place = shuffle_in_place;
+
+ ss.sample = sample;
+
+ ss.sample_covariance = sample_covariance;
+ ss.sample_correlation = sample_correlation;
+ ss.sample_variance = sample_variance;
+ ss.sample_standard_deviation = sample_standard_deviation;
+ ss.sample_skewness = sample_skewness;
+
+ ss.geometric_mean = geometric_mean;
+ ss.harmonic_mean = harmonic_mean;
+ ss.variance = variance;
+ ss.t_test = t_test;
+ ss.t_test_two_sample = t_test_two_sample;jenks
+ + ss.jenksMatrices = jenksMatrices;
+ ss.jenksBreaks = jenksBreaks;
+ ss.jenks = jenks;
+
+ ss.bayesian = bayesian;Distribution-related methods
+ + ss.epsilon = epsilon; // We make ε available to the test suite.
+ ss.factorial = factorial;
+ ss.bernoulli_distribution = bernoulli_distribution;
+ ss.binomial_distribution = binomial_distribution;
+ ss.poisson_distribution = poisson_distribution;
+ ss.chi_squared_goodness_of_fit = chi_squared_goodness_of_fit;Normal distribution
+ + ss.z_score = z_score;
+ ss.cumulative_std_normal_probability = cumulative_std_normal_probability;
+ ss.standard_normal_table = standard_normal_table;Alias this into its common name
+ + ss.average = mean;
+ ss.interquartile_range = iqr;
+ ss.mixin = mixin;
+ ss.median_absolute_deviation = mad;
+
+})(this);