Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 | 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 201125x 201125x 201125x 183535x 183535x 17590x 201125x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 804x 804x 804x 804x 804x 804x 804x 804x 804x 3x 3x 804x 1x 1x 804x 100x 100x 100x 700x 700x 700x 700x 700x 700x 804x 111099x 111099x 111099x 111099x 111099x 700x 700x 111099x 700x 700x 804x 90026x 90026x 90026x 90026x 90026x 199x 199x 90026x 700x 804x 1x 1x 1x 1x 1x | /**
* @license Apache-2.0
*
* Copyright (c) 2018 The Stdlib Authors.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
'use strict';
// MODULES //
var isnan = require( '@stdlib/math/base/assert/is-nan' );
var floor = require( '@stdlib/math/base/special/floor' );
var log1p = require( '@stdlib/math/base/special/log1p' );
var exp = require( '@stdlib/math/base/special/exp' );
var ln = require( '@stdlib/math/base/special/ln' );
var LN_TWO_PI = require( '@stdlib/constants/float64/ln-two-pi' );
var EPS = require( '@stdlib/constants/float64/eps' );
// VARIABLES //
// Value of `lambda` above which the asymptotic expansion is accurate to within double-precision (relative error < 1e-16):
var THRESHOLD = 5000.0;
// Relative weight below which the probability mass of a term no longer affects the sum:
var TOL = EPS * 1.0e-3;
// FUNCTIONS //
/**
* Returns `ln(a/b)` for positive `a` and `b`.
*
* ## Notes
*
* - When `a` and `b` are within a factor of two of each other, `a-b` is exact (Sterbenz lemma), and `log1p` avoids the cancellation which `ln(a/b)` suffers when `a/b` is close to one.
*
* @private
* @param {PositiveNumber} a - numerator
* @param {PositiveNumber} b - denominator
* @returns {number} natural logarithm of the ratio
*/
function lnratio( a, b ) {
var r = a / b;
if ( r > 0.5 && r < 2.0 ) {
return log1p( ( a-b ) / b );
}
return ln( r );
}
// MAIN //
/**
* Returns the entropy of a Poisson distribution.
*
* ## Method
*
* For \\(\lambda < 5000\\), the function evaluates the definition \\(H = -\sum_k p_k \ln p_k\\) using weights relative to the mode \\(m = \lfloor \lambda \rfloor\\), \\(w_k = p_k / p_m\\). As the probabilities sum to one, \\(p_m = 1/W\\), where \\(W = \sum_k w_k\\), and thus
*
* ```tex
* H = \ln W - \frac{1}{W} \sum_k w_k \ln w_k
* ```
*
* The logarithms \\(\ln w_k\\) are accumulated using the recurrences \\(w_{k+1} = w_k \lambda/(k+1)\\) and \\(w_{k-1} = w_k k/\lambda\\). This avoids evaluating \\(\lambda^k/k!\\), which overflows for large \\(k\\), and avoids the cancellation between \\(k \ln \lambda\\), \\(\lambda\\), and \\(\ln k!\\), each of which is much larger than the entropy.
*
* For \\(\lambda \geq 5000\\), the function uses the asymptotic expansion
*
* ```tex
* H = \frac{1}{2} \ln(2 \pi e \lambda) - \frac{1}{12\lambda} - \frac{1}{24\lambda^2} - \frac{19}{360\lambda^3} + O\left(\lambda^{-4}\right)
* ```
*
* whose truncation error is less than \\(10^{-16}\\) relative to the result.
*
* @param {NonNegativeNumber} lambda - mean parameter
* @returns {PositiveNumber} entropy
*
* @example
* var v = entropy( 9.0 );
* // returns ~2.508
*
* @example
* var v = entropy( 1.0 );
* // returns ~1.305
*
* @example
* var v = entropy( -0.2 );
* // returns NaN
*
* @example
* var v = entropy( NaN );
* // returns NaN
*/
function entropy( lambda ) {
var sum;
var wm1;
var lw;
var l2;
var w;
var m;
var k;
if ( isnan( lambda ) || lambda < 0.0 ) {
return NaN;
}
if ( lambda === 0.0 ) {
return 0.0;
}
if ( lambda >= THRESHOLD ) {
l2 = lambda * lambda;
return ( 0.5 * ( LN_TWO_PI + 1.0 + ln( lambda ) ) ) - ( 1.0 / ( 12.0*lambda ) ) - ( 1.0 / ( 24.0*l2 ) ) - ( 19.0 / ( 360.0*l2*lambda ) ); // eslint-disable-line max-len
}
m = floor( lambda );
wm1 = 0.0; // W - 1
sum = 0.0; // Σ w_k ln(w_k)
// Terms above the mode, using w_{k+1} = w_k * λ/(k+1):
lw = 0.0;
for ( k = m; ; k++ ) {
lw += lnratio( lambda, k+1 );
w = exp( lw );
wm1 += w;
sum += w * lw;
if ( w < TOL ) {
break;
}
}
// Terms below the mode, using w_{k-1} = w_k * k/λ:
lw = 0.0;
for ( k = m; k > 0; k-- ) {
lw += lnratio( k, lambda );
w = exp( lw );
wm1 += w;
sum += w * lw;
if ( w < TOL ) {
break;
}
}
return log1p( wm1 ) - ( sum / ( 1.0+wm1 ) );
}
// EXPORTS //
module.exports = entropy;
|