Skip to content

Commit 641243f

Browse files
committed
update specs
1 parent 25bf66f commit 641243f

File tree

2 files changed

+154
-0
lines changed

2 files changed

+154
-0
lines changed

doc/specs/stdlib_stats.md

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -283,6 +283,110 @@ If `mask` is specified, the result is the _k_-th (central) moment of all elemen
283283
{!example/stats/example_moment.f90!}
284284
```
285285

286+
## `pca` - Principal Component Analysis
287+
288+
### Status
289+
290+
Experimental
291+
292+
### Description
293+
294+
Performs Principal Component Analysis (PCA) on a 2D array of observations and features.
295+
The input matrix `x` has shape `(m, n)`, where `m` is the number of observations and `n` is the number of features.
296+
The subroutine computes the principal components (loadings), the singular values, and optionally the feature means.
297+
298+
Two methods are supported:
299+
- `"svd"`: (Default) Computes PCA via Singular Value Decomposition of the centered data. This is generally more numerically stable.
300+
- `"eig"` or `"cov"`: Computes PCA via Eigendecomposition of the covariance matrix.
301+
302+
### Syntax
303+
304+
`call ` [[stdlib_stats(module):pca(interface)]] `(x, components, singular_values [, x_mean [, method [, overwrite_x [, err]]]])`
305+
306+
### Class
307+
308+
Generic subroutine
309+
310+
### Arguments
311+
312+
`x`: Shall be a rank-2 real array with shape `(m, n)`. It is an `intent(inout)` argument. If `overwrite_x` is `.true.`, `x` may be modified during computation.
313+
314+
`components`: Shall be a rank-2 real array with shape `(n_components, n)`. It stores the principal components as rows. It is an `intent(out)` argument.
315+
316+
`singular_values`: Shall be a rank-1 real array with shape `(n_components)`. It stores the singular values in descending order. It is an `intent(out)` argument.
317+
318+
`x_mean` (optional): Shall be a rank-1 real array with shape `(n)`. It stores the mean of each feature (column). It is an `intent(out)` argument.
319+
320+
`method` (optional): Shall be a character string. Either `"svd"` or `"eig"`/`"cov"`. It is an `intent(in)` argument.
321+
322+
`overwrite_x` (optional): Shall be a scalar of type `logical`. If `.true.`, the input matrix `x` can be used as a workspace and modified. It is an `intent(in)` argument.
323+
324+
`err` (optional): Shall be of type `linalg_state_type`. It is an `intent(out)` argument.
325+
326+
### Example
327+
328+
```fortran
329+
{!example/stats/example_pca.f90!}
330+
```
331+
332+
## `pca_transform` - Projects data into principal component space
333+
334+
### Status
335+
336+
Experimental
337+
338+
### Description
339+
340+
Projects the input data `x` into the reduced dimensional space defined by the provided principal components.
341+
The transformation is defined as `x_transformed = (x - x_mean) * components^T`.
342+
343+
### Syntax
344+
345+
`call ` [[stdlib_stats(module):pca_transform(interface)]] `(x, components [, x_mean], x_transformed)`
346+
347+
### Class
348+
349+
Generic subroutine
350+
351+
### Arguments
352+
353+
`x`: Shall be a rank-2 real array with shape `(m, n)`. It is an `intent(in)` argument.
354+
355+
`components`: Shall be a rank-2 real array with shape `(nc, n)`. It stores the principal components as rows. It is an `intent(in)` argument.
356+
357+
`x_mean` (optional): Shall be a rank-1 real array with shape `(n)`. It stores the feature means to subtract. It is an `intent(in)` argument.
358+
359+
`x_transformed`: Shall be a rank-2 real array with shape `(m, nc)`. It stores the projected data. It is an `intent(out)` argument.
360+
361+
## `pca_inverse_transform` - Reconstructs original data from principal component space
362+
363+
### Status
364+
365+
Experimental
366+
367+
### Description
368+
369+
Reconstructs the original data representation from the principal component space.
370+
The reconstruction is defined as `x_reconstructed = x_reduced * components + x_mean`.
371+
372+
### Syntax
373+
374+
`call ` [[stdlib_stats(module):pca_inverse_transform(interface)]] `(x_reduced, components [, x_mean], x_reconstructed)`
375+
376+
### Class
377+
378+
Generic subroutine
379+
380+
### Arguments
381+
382+
`x_reduced`: Shall be a rank-2 real array with shape `(m, nc)`. It is an `intent(in)` argument.
383+
384+
`components`: Shall be a rank-2 real array with shape `(nc, n)`. It stores the principal components as rows. It is an `intent(in)` argument.
385+
386+
`x_mean` (optional): Shall be a rank-1 real array with shape `(n)`. It stores the feature means to add back. It is an `intent(in)` argument.
387+
388+
`x_reconstructed`: Shall be a rank-2 real array with shape `(m, n)`. It stores the reconstructed data. It is an `intent(out)` argument.
389+
286390
## `var` - variance of array elements
287391

288392
### Status

example/stats/example_pca.f90

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
program example_pca
2+
use stdlib_kinds, only: dp
3+
use stdlib_stats, only: pca, pca_transform, pca_inverse_transform
4+
use stdlib_linalg_state, only: linalg_state_type
5+
implicit none
6+
7+
real(dp) :: x(3, 2), components(2, 2), s(2), mu(2)
8+
real(dp) :: x_trans(3, 2), x_inv(3, 2)
9+
type(linalg_state_type) :: err
10+
integer :: i
11+
12+
! Input data: 3 observations, 2 features
13+
x = reshape([1.0_dp, 3.0_dp, 5.0_dp, 2.0_dp, 4.0_dp, 6.0_dp], [3, 2])
14+
15+
print *, "Original data:"
16+
do i = 1, 3
17+
print "(2f6.2)", x(i, :)
18+
end do
19+
20+
! Perform PCA
21+
call pca(x, components, s, x_mean=mu, err=err)
22+
23+
if (err%ok()) then
24+
print *, ""
25+
print *, "Feature means:", mu
26+
print *, "Singular values:", s
27+
print *, "Principal components (as rows):"
28+
print "(2f6.3)", components(1, :)
29+
print "(2f6.3)", components(2, :)
30+
31+
! Transform data to principal components space
32+
call pca_transform(x, components, mu, x_trans)
33+
print *, ""
34+
print *, "Transformed data (projected):"
35+
do i = 1, 3
36+
print "(2f8.3)", x_trans(i, :)
37+
end do
38+
39+
! Inverse transform to reconstruct original data
40+
call pca_inverse_transform(x_trans, components, mu, x_inv)
41+
print *, ""
42+
print *, "Reconstructed data:"
43+
do i = 1, 3
44+
print "(2f6.2)", x_inv(i, :)
45+
end do
46+
else
47+
print *, "PCA failed: ", err%message
48+
end if
49+
50+
end program example_pca

0 commit comments

Comments
 (0)