|
16 | 16 | * limitations under the License.
|
17 | 17 | */
|
18 | 18 |
|
19 |
| -#include "stdlib/stats/base/dists/triangular/mgf.h" |
20 | 19 | #include "stdlib/math/base/assert/is_nan.h"
|
21 | 20 | #include "stdlib/math/base/special/exp.h"
|
22 |
| -#include "stdlib/math/base/special/pow.h" |
| 21 | +#include "stdlib/math/base/special/expm1.h" |
| 22 | +#include "stdlib/stats/base/dists/triangular/mgf.h" |
| 23 | + |
| 24 | +/** |
| 25 | +* Helper function for repeated computation in the MGF formula. |
| 26 | +* |
| 27 | +* @private |
| 28 | +* @param x input value |
| 29 | +* @return evaluated result |
| 30 | +*/ |
| 31 | +static double phi2( const double x ) { |
| 32 | + if ( x == 0.0 ) { |
| 33 | + return 1.0; |
| 34 | + } |
| 35 | + return ( 2.0 * ( stdlib_base_expm1( x ) - x ) ) / ( x * x ); |
| 36 | +} |
23 | 37 |
|
24 | 38 | /**
|
25 | 39 | * Evaluates the moment-generating function (MGF) for a triangular distribution with lower limit `a`, upper limit `b`, and mode `c` at a value `t`.
|
|
35 | 49 | * // returns ~1.021
|
36 | 50 | */
|
37 | 51 | double stdlib_base_dists_triangular_mgf( const double t, const double a, const double b, const double c ) {
|
38 |
| - double bmc; |
39 |
| - double bma; |
40 |
| - double cma; |
41 |
| - double ret; |
42 |
| - if ( |
43 |
| - stdlib_base_is_nan( t ) || |
44 |
| - stdlib_base_is_nan( a ) || |
45 |
| - stdlib_base_is_nan( b ) || |
46 |
| - stdlib_base_is_nan( c ) || |
47 |
| - a > c || |
48 |
| - c > b |
49 |
| - ) { |
50 |
| - return 0.0/0.0; // NaN |
| 52 | + if ( stdlib_base_is_nan( t ) || stdlib_base_is_nan( a ) || stdlib_base_is_nan( b ) || stdlib_base_is_nan( c ) || a > c || c > b ) { |
| 53 | + return 0.0 / 0.0; // NaN |
51 | 54 | }
|
52 | 55 | if ( t == 0.0 ) {
|
53 | 56 | return 1.0;
|
54 | 57 | }
|
55 |
| - bmc = b - c; |
56 |
| - bma = b - a; |
57 |
| - cma = c - a; |
58 |
| - ret = ( bmc * stdlib_base_exp( a*t ) ) - ( bma * stdlib_base_exp( c*t ) ); |
59 |
| - ret += cma * stdlib_base_exp( b*t ); |
60 |
| - ret *= 2.0; |
61 |
| - ret /= bma * cma * bmc * stdlib_base_pow( t, 2.0 ); |
62 |
| - return ret; |
| 58 | + |
| 59 | + double term1 = 0.0; |
| 60 | + double term2 = 0.0; |
| 61 | + |
| 62 | + if ( a < c ) { |
| 63 | + if ( c < b ) { |
| 64 | + term1 = ( c - a ) * phi2( ( a - c ) * t ); |
| 65 | + term2 = ( b - c ) * phi2( ( b - c ) * t ); |
| 66 | + return stdlib_base_exp( c * t ) * ( term1 + term2 ) / ( b - a ); |
| 67 | + } |
| 68 | + return stdlib_base_exp( c * t ) * phi2( ( a - c ) * t ); |
| 69 | + } |
| 70 | + |
| 71 | + if ( c < b ) { |
| 72 | + return stdlib_base_exp( c * t ) * phi2( ( b - c ) * t ); |
| 73 | + } |
| 74 | + |
| 75 | + return stdlib_base_exp( c * t ); |
63 | 76 | }
|
0 commit comments