SF_STAT

Description

SciFortran module for statistical estimators and visualization

Quick access

Types:

histogram, pdf_kernel, pdf_kernel_2d

Routines:

get_covariance(), get_curt(), get_mean(), get_moments(), get_sd(), get_skew(), get_var(), histogram_accumulate(), histogram_allocate(), histogram_deallocate(), histogram_get_range(), histogram_get_value(), histogram_print(), histogram_reset(), histogram_set_range_uniform(), pdf_accumulate(), pdf_allocate(), pdf_curt(), pdf_deallocate(), pdf_mean(), pdf_moment(), pdf_normalize(), pdf_print(), pdf_print_moments(), pdf_push_sigma(), pdf_read(), pdf_save(), pdf_sdev(), pdf_set_range(), pdf_sigma(), pdf_skew(), pdf_var()

Used modules

Types

type  sf_stat/histogram

This derived type represents a histogram of n bins. The bin \(i\), with \(i=0,\dots,n-1\), covers the interval [range(i),range(i+1)) and contains in bin(i) the sum of the weights of the values accumulated in it. It is created by histogram_allocate and handled by the histogram_* routines.

Type fields:

  • n [integer] – number of bins

  • range (•) [real, pointer] – edges of the bins, range(0:n); bin i is [range(i),range(i+1))

  • bin (•) [real, pointer] – content of the bins, bin(0:n); only bin(0:n-1) is used

type  sf_stat/pdf_kernel

This derived type holds the probability density function (pdf) of a one-dimensional variable, estimated with the Gaussian kernel density method. The pdf is represented on a uniform grid x of N points between xmin and xmax: every accumulated data point adds to it a Gaussian of width sigma centered on the data point. See pdf_allocate, pdf_set_range, pdf_accumulate and pdf_normalize.

Type fields:

  • n [integer] – number of grid points

  • xmin [real] – lower bound of the grid

  • xmax [real] – upper bound of the grid

  • dx [real] – grid step

  • ndata [integer] – number of data points accumulated

  • x (•) [real, allocatable] – grid, size N

  • pdf (•) [real, allocatable] – pdf on the grid, size N

  • sigma [real] – width (standard deviation) of the Gaussian kernel

  • status [logical] – .true. if the object is allocated

  • variance [logical] – .true. if sigma has been stored with pdf_push_sigma

  • rescale [logical] – .true. once the pdf has been normalized

type  sf_stat/pdf_kernel_2d

This derived type holds the probability density function (pdf) of a two-dimensional variable, estimated with the Gaussian kernel density method. The pdf is represented on a uniform grid of N(1) x N(2) points, defined by the grids x and y between xmin and xmax: every accumulated data point adds to it a bivariate Gaussian with covariance matrix Sigma centered on the data point. See pdf_allocate, pdf_set_range, pdf_accumulate and pdf_normalize.

Type fields:

  • n (2) [integer] – number of grid points along x and y

  • xmin (2) [real] – lower bounds of the grid, (x,y)

  • xmax (2) [real] – upper bounds of the grid, (x,y)

  • dx (2) [real] – grid steps, (x,y)

  • ndata [integer] – number of data points accumulated

  • x (•) [real, allocatable] – grid along x, size N(1)

  • y (•) [real, allocatable] – grid along y, size N(2)

  • pdf (•, •) [real, allocatable] – pdf on the grid, [N(1),N(2)]

  • sigma (2, 2) [real] – covariance matrix of the Gaussian kernel

  • status [logical] – .true. if the object is allocated

  • variance [logical] – .true. if Sigma has been stored with pdf_push_sigma

  • rescale [logical] – .true. once the pdf has been normalized

Subroutines and functions

interface  sf_stat/pdf_allocate(self, n, nvec)

Generic interface to allocate a pdf object: a pdf_kernel with N grid points (pdf_allocate_1d), or a pdf_kernel_2d with a grid of Nvec(1) x Nvec(2) points (pdf_allocate_2d). The pdf and the kernel width are set to zero, no data is accumulated, and the object is flagged as allocated. The grid is not defined until pdf_set_range is called. The program stops if the object is already allocated.

Parameters:

  • self [pdf_kernel] – pdf object, not allocated

  • n [integer] – number of grid points

  • nvec (2) [integer] – number of grid points along x and y

interface  sf_stat/pdf_deallocate(self)

Generic interface to deallocate a pdf object, pdf_deallocate_1d or pdf_deallocate_2d: the arrays are freed and the ranges, the kernel width and the number of data are reset to zero. The program stops if the object is not allocated.

Parameters:

self [pdf_kernel] – pdf object

interface  sf_stat/pdf_save(self, pfile)

Generic interface to save a pdf object, pdf_save_1d or pdf_save_2d, to the formatted file pfile. The file is overwritten, and it contains one record per item: the number of grid points, the ranges, the grid steps, the number of data, the grids, the pdf, the kernel width and the three flags of the object. It can be read back with pdf_read. The program stops if the object is not allocated.

Parameters:

  • self [pdf_kernel] – pdf object

  • pfile [character(len=*)] – name of the output file

interface  sf_stat/pdf_read(self, pfile)

Generic interface to read a pdf object, pdf_read_1d or pdf_read_2d, from a file written by pdf_save. The object is allocated by the routine, therefore it must not be allocated on input, otherwise the program stops.

Parameters:

  • self [pdf_kernel] – pdf object, not allocated; allocated and set on output

  • pfile [character(len=*)] – name of the file written by pdf_save

interface  sf_stat/pdf_set_range(self, a, b)

Generic interface to define the grid of a pdf object, as a uniform mesh between the lower bound a and the upper bound b: a scalar pair for pdf_set_range_1d, vectors with the bounds along x and y for pdf_set_range_2d. It sets the ranges and the grid steps of the object. The program stops if the object is not allocated.

Parameters:

  • self [pdf_kernel] – pdf object

  • a (various shapes) [real] – lower bound of the grid

  • b (various shapes) [real] – upper bound of the grid

interface  sf_stat/pdf_push_sigma(self, sigma)

Generic interface to store the width of the Gaussian kernel in the pdf object: the standard deviation sigma in one dimension (pdf_push_sigma_1d), the \(2\times 2\) covariance matrix in two dimensions (pdf_push_sigma_2d). Once stored, it is used by pdf_accumulate, and the argument sigma of the latter is ignored. The program stops if the object is not allocated.

Parameters:

  • self [pdf_kernel] – pdf object

  • sigma (various shapes) [real] – width (standard deviation) of the Gaussian kernel

interface  sf_stat/pdf_sigma(self, data, h, sdev, n, nvec)

Generic interface to estimate the width of the Gaussian kernel with Silverman’s rule of thumb. The result is returned in h, the pdf object is not modified (but it must be allocated). In one dimension \(h = (4/3N)^{1/5}\sigma\), with \(N\) the number of data and \(\sigma\) their standard deviation; in two dimensions the result is the diagonal covariance matrix \(h_{ii} = N^{-1/3}\sigma_i^2\).

The standard deviation is computed from the data (pdf_sigma_data_1d, and pdf_sigma_data_2d for data of shape [2,L]) or it is given together with the number of data (pdf_sigma_sdev_1d, pdf_sigma_sdev_2d). The program stops if the object is not allocated or, in two dimensions, if the data do not have shape [2,L].

Parameters:

  • self [pdf_kernel] – pdf object

  • data (various shapes) [real] – sample of data

  • h (various shapes) [real] – width of the kernel, Silverman rule

  • sdev (various shapes) [real] – standard deviation of the data

  • n [integer] – number of data

  • nvec (2) [integer] – number of data, (x,y)

interface  sf_stat/pdf_accumulate(self, data[, sigma])

Generic interface to accumulate data in a pdf object. Every data point adds to the pdf on the grid a Gaussian kernel centered on it, and the number of data of the object is increased. In one dimension the kernel is \(e^{-(x-m)^2/2\sigma^2}/\sqrt{2\pi}\sigma\) and the data is a scalar (pdf_accumulate_s_1d) or a vector, whose elements are accumulated one by one (pdf_accumulate_v_1d). In two dimensions it is \(e^{-\frac{1}{2}(x-m)^T\Sigma^{-1}(x-m)}/2\pi\sqrt{\det\Sigma}\) and the data is a point (x,y) (pdf_accumulate_s_2d).

The width of the kernel is the one stored with pdf_push_sigma; if none has been stored it must be passed in the optional argument sigma, otherwise the program stops. The program stops also if the object is not allocated.

Parameters:

  • self [pdf_kernel] – pdf object

  • data (various shapes) [real] – data point

Options:

sigma (various shapes) [real] – optional: kernel width, if not stored in self

interface  sf_stat/pdf_normalize(self)

Generic interface to normalize a pdf object, pdf_normalize_1d or pdf_normalize_2d. The pdf is divided by the number of accumulated data (the first time only) and then rescaled so that its integral on the grid is one, i.e. \(\sum_i p_i \Delta x = 1\) (\(\sum_{ij} p_{ij} \Delta x \Delta y = 1\) in two dimensions). It is called by pdf_print unless otherwise requested. The program stops if the object is not allocated.

Parameters:

self [pdf_kernel] – pdf object

interface  sf_stat/pdf_print(self, pfile[, normalize])

Generic interface to write a pdf object to the file pfile. The pdf is normalized first, with pdf_normalize, unless normalize is .false.. In one dimension the two columns x pdf are appended to the file, followed by a blank line (pdf_print_pfile_1d); in two dimensions the pdf is written with splot3d (pdf_print_pfile_2d). The program stops if the object is not allocated.

Parameters:

  • self [pdf_kernel] – pdf object

  • pfile [character(len=*)] – name of the output file (appended)

Options:

normalize [logical] – optional: normalize first, default .true.

interface  sf_stat/pdf_mean(self)

Generic interface for the mean of a one-dimensional pdf, \(\int x\,p(x)\,dx\), integrated with the Simpson’s rule on the grid (pdf_mean_1d).

Parameters:

self [pdf_kernel] – pdf object

Result:

mean [real] – mean of the pdf

interface  sf_stat/pdf_var(self)

Generic interface for the variance of a one-dimensional pdf, \(\int (x-\mu)^2\,p(x)\,dx\), integrated with the Simpson’s rule on the grid (pdf_var_1d).

Parameters:

self [pdf_kernel] – pdf object

Result:

var [real] – variance of the pdf

interface  sf_stat/pdf_sdev(self)

Generic interface for the standard deviation of a one-dimensional pdf, the square root of pdf_var (pdf_sdev_1d).

Parameters:

self [pdf_kernel] – pdf object

Result:

sdev [real] – standard deviation of the pdf

interface  sf_stat/pdf_moment(self, n[, mu])

Generic interface for the moment of order \(n\) of a one-dimensional pdf around a point \(\mu\), \(\int (x-\mu)^n\,p(x)\,dx\), integrated with the Simpson’s rule on the grid (pdf_moment_1d). The point \(\mu\) is optional and it is zero by default, i.e. the moment is the raw one.

Parameters:

  • self [pdf_kernel] – pdf object

  • n [integer] – order of the moment

Options:

mu [real] – optional: center of the moment, default 0

Result:

mom [real] – moment of order n of the pdf

interface  sf_stat/pdf_skew(self)

Generic interface for the skewness of a one-dimensional pdf, the third standardized moment \(\int (x-\mu)^3 p(x)\,dx/\sigma^3\) (pdf_skew_1d).

Parameters:

self [pdf_kernel] – pdf object

Result:

skew [real] – skewness of the pdf

interface  sf_stat/pdf_curt(self)

Generic interface for the kurtosis of a one-dimensional pdf, the fourth standardized moment \(\int (x-\mu)^4 p(x)\,dx/\sigma^4\), NOT the excess kurtosis (pdf_curt_1d).

Parameters:

self [pdf_kernel] – pdf object

Result:

curt [real] – kurtosis of the pdf (not excess)

interface  sf_stat/pdf_print_moments(self, pfile)

Generic interface to write to the file pfile the moments of a one-dimensional pdf, on a single line: mean, standard deviation, skewness and kurtosis (pdf_print_moments_pfile_1d).

Parameters:

  • self [pdf_kernel] – pdf object

  • pfile [character(len=*)] – name of the output file

subroutine  sf_stat/get_moments(data[, ave, sdev, var, skew, curt])

This subroutine computes the mean, the standard deviation, the variance, the skewness and the kurtosis of the real array data. All the outputs are optional: only the requested ones are returned.

The variance is the unbiased estimator, with \(n-1\) in the denominator, computed with the corrected two-pass algorithm, and the standard deviation is its square root. The skewness is the third central moment divided by \(n\sigma^3\). The kurtosis is the EXCESS kurtosis, the fourth central moment divided by \(n\sigma^4\) minus 3, which is zero for a normal distribution. If the variance is zero the skewness and the kurtosis are set to zero. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Options:

  • ave [real, out] – mean

  • sdev [real, out] – standard deviation (unbiased)

  • var [real, out] – variance (unbiased)

  • skew [real, out] – skewness

  • curt [real, out] – excess kurtosis

function  sf_stat/get_mean(data)

This function returns the mean of the real array data, computed with get_moments. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Result:

mean [real] – mean

function  sf_stat/get_sd(data)

This function returns the standard deviation, with \(n-1\) in the denominator, of the real array data, computed with get_moments. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Result:

sd [real] – standard deviation

function  sf_stat/get_var(data)

This function returns the variance, with \(n-1\) in the denominator, of the real array data, computed with get_moments. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Result:

var [real] – variance

function  sf_stat/get_skew(data)

This function returns the skewness of the real array data, computed with get_moments. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Result:

skew [real] – skewness

function  sf_stat/get_curt(data)

This function returns the excess kurtosis of the real array data, computed with get_moments. The program stops if data has less than two elements.

Parameters:

data (•) [real, in] – sample, size n >= 2

Result:

curt [real] – excess kurtosis

function  sf_stat/get_covariance(data, mean)

This function returns the covariance matrix of the variables stored as the rows of data, for \(L\) samples stored as the columns: \(C_{ij} = \frac{1}{L-1}\sum_k (x_{ik}-\bar{x}_i)(x_{jk}-\bar{x}_j)\), with \(\bar{x}_i\) the means given in mean.

Parameters:

  • data (•, •) [real, in,required] – samples, [ncol,L]: ncol variables (rows), L samples (columns)

  • mean (size(data, 1)) [real, in] – mean of each variable, size ncol

Result:

covariance (size(data, 1), size(data, 1)) [real] – covariance matrix, [ncol,ncol]

function  sf_stat/histogram_allocate(n)

This function allocates a histogram of n bins, with all the edges and the bin contents set to zero. The edges can be defined with histogram_set_range_uniform. The program stops if n is not positive.

Parameters:

n [integer, in] – number of bins

Result:

h [histogram] – histogram, allocated

subroutine  sf_stat/histogram_deallocate(h)

This subroutine deallocates the histogram h and sets its number of bins to zero.

Parameters:

h [histogram] – histogram to be deallocated

subroutine  sf_stat/histogram_reset(h)

This subroutine sets to zero both the edges and the contents of the bins of the histogram h, which is kept allocated. The edges must therefore be defined again with histogram_set_range_uniform.

Parameters:

h [histogram] – histogram to be reset

subroutine  sf_stat/histogram_set_range_uniform(h, xmin, xmax)

This subroutine defines the edges of the bins of the histogram h as n+1 equally spaced points between xmin and xmax, and it resets the contents of the bins to zero. The program stops if xmin>=xmax.

Parameters:

  • h [histogram, inout] – histogram, allocated

  • xmin [real, in] – lower edge of the first bin

  • xmax [real, in] – upper edge of the last bin

subroutine  sf_stat/histogram_accumulate(h, x, w)

This subroutine adds the weight w to the bin of the histogram h containing the value x, found with a bisection search. If x is outside the range of the histogram the message X out of range! is printed and, as the routine is written, the weight is added to the first bin.

Parameters:

  • h [histogram, inout] – histogram, with the edges defined

  • x [real, in] – value to be accumulated

  • w [real, in] – weight of the value

subroutine  sf_stat/histogram_get_range(h, index, lower, upper)

This subroutine returns the edges of the bin number index of the histogram h: the bins are numbered from 0 to n-1. The program stops if index is not smaller than n.

Parameters:

  • h [histogram, in] – histogram

  • index [integer, in] – number of the bin, from 0 to n-1

  • lower [real, out] – lower edge of the bin

  • upper [real, out] – upper edge of the bin

function  sf_stat/histogram_get_value(h, index)

This function returns the content of the bin number index of the histogram h: the bins are numbered from 0 to n-1. The program stops if index is not smaller than n.

Parameters:

  • h [histogram, in] – histogram

  • index [integer, in] – number of the bin, from 0 to n-1

Result:

value [real] – content of the bin

subroutine  sf_stat/histogram_print(h, unit)

This subroutine writes the histogram h on the unit unit, as a step profile ready to be plotted: pairs edge content with the two edges of each bin, preceded and followed by a point of zero content, in the format 2F12.7.

Parameters:

  • h [histogram, in] – histogram

  • unit [integer, in] – output unit, already open