Skip to content

Commit 9797449

Browse files
committed
feat: add support for skew normal distribution PDF
--- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: passed - task: lint_package_json status: passed - task: lint_repl_help status: passed - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: passed - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: passed - task: lint_python status: na - task: lint_r status: passed - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 59f0fbb commit 9797449

File tree

17 files changed

+1317
-0
lines changed

17 files changed

+1317
-0
lines changed

etc/r/requirements.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,3 +3,4 @@ microbenchmark
33
jsonlite
44
VGAM
55
distr
6+
sn
Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
2+
<!--
3+
4+
@license Apache-2.0
5+
6+
Copyright (c) 2025 The Stdlib Authors.
7+
8+
Licensed under the Apache License, Version 2.0 (the "License");
9+
you may not use this file except in compliance with the License.
10+
You may obtain a copy of the License at
11+
12+
http://www.apache.org/licenses/LICENSE-2.0
13+
14+
Unless required by applicable law or agreed to in writing, software
15+
distributed under the License is distributed on an "AS IS" BASIS,
16+
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
See the License for the specific language governing permissions and
18+
limitations under the License.
19+
20+
-->
21+
22+
# Probability Density Function
23+
24+
> [Skew normal][skew-normal-distribution] distribution probability density function (PDF).
25+
26+
<section class="intro">
27+
28+
The [probability density function][pdf] (PDF) for a [skew normal][skew-normal-distribution] random variable is
29+
30+
```math
31+
f(x;\mu,\sigma,\alpha)=\frac{2}{\sigma}\phi \left({\frac {x-\mu }{\sigma }}\right)\Phi \left(\alpha \left({\frac {x-\mu }{\sigma }}\right)\right)
32+
```
33+
34+
where,
35+
36+
> `ϕ` is the standard normal probability distribution function,
37+
> `Φ` is the standard normal cumulative distribution function,
38+
> `μ` is the mean,
39+
> `σ` is the standard deviation,
40+
> `α` is the skewness.
41+
42+
</section>
43+
44+
<!-- /.intro -->
45+
46+
<section class="usage">
47+
48+
## Usage
49+
50+
```javascript
51+
var pdf = require( '@stdlib/stats/base/dists/skew-normal/pdf' );
52+
```
53+
54+
#### pdf( x, mu, sigma, alpha )
55+
56+
Evaluates the [probability density function][pdf] (PDF) for a [skew normal][skew-normal-distribution] distribution with parameters `mu` (mean), `sigma` (standard deviation), and `alpha` (skewness).
57+
58+
```javascript
59+
var y = pdf( 2.0, 0.0, 1.0, 2.0 );
60+
// returns ~0.108
61+
62+
y = pdf( -1.0, 4.0, 2.0, 0.0 );
63+
// returns ~0.009
64+
```
65+
66+
If provided `NaN` as any argument, the function returns `NaN`.
67+
68+
```javascript
69+
var y = pdf( NaN, 0.0, 1.0, 0.0 );
70+
// returns NaN
71+
72+
y = pdf( 0.0, NaN, 1.0, 0.0 );
73+
// returns NaN
74+
75+
y = pdf( 0.0, 0.0, NaN, 0.0 );
76+
// returns NaN
77+
78+
y = pdf( 0.0, 0.0, 0.0, NaN );
79+
// returns NaN
80+
```
81+
82+
If provided `sigma < 0`, the function returns `NaN`.
83+
84+
```javascript
85+
var y = pdf( 2.0, 0.0, -1.0, 0.0 );
86+
// returns NaN
87+
```
88+
89+
If provided `sigma = 0`, the function evaluates the [PDF][pdf] of a [degenerate distribution][degenerate-distribution] centered at `mu`.
90+
91+
```javascript
92+
var y = pdf( 2.0, 8.0, 0.0, 0.0 );
93+
// returns 0.0
94+
95+
y = pdf( 8.0, 8.0, 0.0, 0.0 );
96+
// returns Infinity
97+
```
98+
99+
#### pdf.factory( mu, sigma, alpha )
100+
101+
Partially apply `mu`, `sigma`, and `alpha` to create a reusable `function` for evaluating the PDF.
102+
103+
```javascript
104+
var mypdf = pdf.factory( 10.0, 2.0, -1.0 );
105+
106+
var y = mypdf( 10.0 );
107+
// returns ~0.199
108+
109+
y = mypdf( 5.0 );
110+
// returns ~0.017
111+
```
112+
113+
</section>
114+
115+
<!-- /.usage -->
116+
117+
<section class="examples">
118+
119+
## Examples
120+
121+
<!-- eslint no-undef: "error" -->
122+
123+
```javascript
124+
var randu = require( '@stdlib/random/base/randu' );
125+
var pdf = require( '@stdlib/stats/base/dists/skew-normal/pdf' );
126+
127+
var alpha;
128+
var sigma;
129+
var mu;
130+
var x;
131+
var y;
132+
var i;
133+
134+
for ( i = 0; i < 10; i++ ) {
135+
x = randu() * 10.0;
136+
mu = (randu() * 10.0) - 5.0;
137+
sigma = randu() * 20.0;
138+
alpha = (randu() * 40.0) - 20.0;
139+
y = pdf( x, mu, sigma, alpha );
140+
console.log( 'x: %d, µ: %d, σ: %d, α: %d, f(x;µ,σ,α): %d', x, mu, sigma, alpha, y );
141+
}
142+
```
143+
144+
</section>
145+
146+
<!-- /.examples -->
147+
148+
<!-- Section for related `stdlib` packages. Do not manually edit this section, as it is automatically populated. -->
149+
150+
<section class="related">
151+
152+
</section>
153+
154+
<!-- /.related -->
155+
156+
<!-- Section for all links. Make sure to keep an empty line after the `section` element and another before the `/section` close. -->
157+
158+
<section class="links">
159+
160+
[pdf]: https://en.wikipedia.org/wiki/Probability_density_function
161+
162+
[skew-normal-distribution]: https://en.wikipedia.org/wiki/Skew_normal_distribution
163+
164+
[degenerate-distribution]: https://en.wikipedia.org/wiki/Degenerate_distribution
165+
166+
</section>
167+
168+
<!-- /.links -->
Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2025 The Stdlib Authors.
5+
*
6+
* Licensed under the Apache License, Version 2.0 (the "License");
7+
* you may not use this file except in compliance with the License.
8+
* You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing, software
13+
* distributed under the License is distributed on an "AS IS" BASIS,
14+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
* See the License for the specific language governing permissions and
16+
* limitations under the License.
17+
*/
18+
19+
'use strict';
20+
21+
// MODULES //
22+
23+
var bench = require( '@stdlib/bench' );
24+
var randu = require( '@stdlib/random/base/randu' );
25+
var isnan = require( '@stdlib/math/base/assert/is-nan' );
26+
var EPS = require( '@stdlib/constants/float64/eps' );
27+
var pkg = require( './../package.json' ).name;
28+
var pdf = require( './../lib' );
29+
30+
31+
// MAIN //
32+
33+
bench( pkg, function benchmark( b ) {
34+
var alpha;
35+
var sigma;
36+
var mu;
37+
var x;
38+
var y;
39+
var i;
40+
41+
b.tic();
42+
for ( i = 0; i < b.iterations; i++ ) {
43+
x = ( randu()*200.0 ) - 100;
44+
mu = ( randu()*100.0 ) - 50.0;
45+
sigma = ( randu()*20.0 ) + EPS;
46+
alpha = ( randu()*40.0 ) - 20.0;
47+
y = pdf( x, mu, sigma, alpha );
48+
if ( isnan( y ) ) {
49+
b.fail( 'should not return NaN' );
50+
}
51+
}
52+
b.toc();
53+
if ( isnan( y ) ) {
54+
b.fail( 'should not return NaN' );
55+
}
56+
b.pass( 'benchmark finished' );
57+
b.end();
58+
});
59+
60+
bench( pkg+':factory', function benchmark( b ) {
61+
var mypdf;
62+
var alpha;
63+
var sigma;
64+
var mu;
65+
var x;
66+
var y;
67+
var i;
68+
69+
mu = 0.0;
70+
sigma = 1.5;
71+
alpha = 0.5;
72+
mypdf = pdf.factory( mu, sigma, alpha);
73+
74+
b.tic();
75+
for ( i = 0; i < b.iterations; i++ ) {
76+
x = ( randu()*6.0 ) - 3.0;
77+
y = mypdf( x );
78+
if ( isnan( y ) ) {
79+
b.fail( 'should not return NaN' );
80+
}
81+
}
82+
b.toc();
83+
if ( isnan( y ) ) {
84+
b.fail( 'should not return NaN' );
85+
}
86+
b.pass( 'benchmark finished' );
87+
b.end();
88+
});
Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
Package: sn-benchmarks
2+
Title: Benchmarks
3+
Version: 0.0.0
4+
Authors@R: person("stdlib", "js", role = c("aut","cre"))
5+
Description: Benchmarks.
6+
Depends: R (>=3.4.0)
7+
Imports:
8+
microbenchmark
9+
sn
10+
LazyData: true
Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#!/usr/bin/env Rscript
2+
#
3+
# @license Apache-2.0
4+
#
5+
# Copyright (c) 2025 The Stdlib Authors.
6+
#
7+
# Licensed under the Apache License, Version 2.0 (the "License");
8+
# you may not use this file except in compliance with the License.
9+
# You may obtain a copy of the License at
10+
#
11+
# http://www.apache.org/licenses/LICENSE-2.0
12+
#
13+
# Unless required by applicable law or agreed to in writing, software
14+
# distributed under the License is distributed on an "AS IS" BASIS,
15+
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16+
# See the License for the specific language governing permissions and
17+
# limitations under the License.
18+
19+
# Set the precision to 16 digits:
20+
options( digits = 16L );
21+
22+
#' Run benchmarks.
23+
#'
24+
#' @examples
25+
#' main();
26+
main <- function() {
27+
# Define benchmark parameters:
28+
name <- 'dist-skew-normal-pdf';
29+
iterations <- 1000000L;
30+
repeats <- 3L;
31+
32+
#' Print the TAP version.
33+
#'
34+
#' @examples
35+
#' print_version();
36+
print_version <- function() {
37+
cat( 'TAP version 13\n' );
38+
}
39+
40+
#' Print the TAP summary.
41+
#'
42+
#' @param total Total number of tests.
43+
#' @param passing Total number of passing tests.
44+
#'
45+
#' @examples
46+
#' print_summary( 3, 3 );
47+
print_summary <- function( total, passing ) {
48+
cat( '#\n' );
49+
cat( paste0( '1..', total, '\n' ) ); # TAP plan
50+
cat( paste0( '# total ', total, '\n' ) );
51+
cat( paste0( '# pass ', passing, '\n' ) );
52+
cat( '#\n' );
53+
cat( '# ok\n' );
54+
}
55+
56+
#' Print benchmark results.
57+
#'
58+
#' @param iterations Number of iterations.
59+
#' @param elapsed Elapsed time in seconds.
60+
#'
61+
#' @examples
62+
#' print_results( 10000L, 0.131009101868 );
63+
print_results <- function( iterations, elapsed ) {
64+
rate <- iterations / elapsed;
65+
cat( ' ---\n' );
66+
cat( paste0( ' iterations: ', iterations, '\n' ) );
67+
cat( paste0( ' elapsed: ', elapsed, '\n' ) );
68+
cat( paste0( ' rate: ', rate, '\n' ) );
69+
cat( ' ...\n' );
70+
}
71+
72+
#' Run a benchmark.
73+
#'
74+
#' ## Notes
75+
#'
76+
#' * We compute and return a total "elapsed" time, rather than the minimum
77+
#' evaluation time, to match benchmark results in other languages (e.g.,
78+
#' Python).
79+
#'
80+
#'
81+
#' @param iterations Number of Iterations.
82+
#' @return Elapsed time in seconds.
83+
#'
84+
#' @examples
85+
#' elapsed <- benchmark( 10000L );
86+
benchmark <- function( iterations ) {
87+
# Run the benchmarks:
88+
results <- microbenchmark::microbenchmark(sn::dsn(
89+
runif( 1.0, -100.0, 100.0 ),
90+
runif( 1.0, -50.0, 50.0 ),
91+
runif( 1.0, .Machine$double.eps, 20.0 ), # nolint
92+
runif( 1.0, -20.0, 20.0 )
93+
), times = iterations );
94+
95+
# Sum all the raw timing results to get a total "elapsed" time:
96+
elapsed <- sum( results$time ); # nolint
97+
98+
# Convert the elapsed time from nanoseconds to seconds:
99+
elapsed <- elapsed / 1.0e9;
100+
101+
return( elapsed );
102+
}
103+
104+
print_version();
105+
for ( i in 1L:repeats ) {
106+
cat( paste0( '# r::', name, '\n' ) );
107+
elapsed <- benchmark( iterations );
108+
print_results( iterations, elapsed );
109+
cat( paste0( 'ok ', i, ' benchmark finished', '\n' ) );
110+
}
111+
print_summary( repeats, repeats );
112+
}
113+
114+
main();

0 commit comments

Comments
 (0)