Package 'ica'

Title: Independent Component Analysis
Description: Independent Component Analysis (ICA) using various algorithms: FastICA, Information-Maximization (Infomax), and Joint Approximate Diagonalization of Eigenmatrices (JADE).
Authors: Nathaniel E. Helwig <[email protected]>
Maintainer: Nathaniel E. Helwig <[email protected]>
License: GPL (>= 2)
Version: 1.0-3
Built: 2024-12-19 06:43:40 UTC
Source: CRAN

Help Index


Amari-Cichocki-Yang Error

Description

The Amari-Cichocki-Yang (ACY) error is an asymmetric measure of dissimilarity between two nonsingular matrices X and Y. The ACY error: (a) is invariant to permutation and rescaling of the columns of X and Y, (b) ranges between 0 and n-1, and (c) equals 0 if and only if X and Y are identical up to column permutations and rescalings.

Usage

acy(X,Y)

Arguments

X

Nonsingular matrix of dimension n×nn \times n (test matrix).

Y

Nonsingular matrix of dimension n×nn \times n (target matrix).

Details

The ACY error is defined as

12ni=1n(j=1naijmaxjaij1)+12nj=1n(i=1naijmaxiaij1)\frac{1}{2n}\sum_{i=1}^{n}\left(\frac{\sum_{j=1}^{n}|a_{ij}|}{\max_{j}|a_{ij}|}-1\right) + \frac{1}{2n}\sum_{j=1}^{n}\left(\frac{\sum_{i=1}^{n}|a_{ij}|}{\max_{i}|a_{ij}|}-1\right)

where aij=(Y1X)ija_{ij} = (\mathbf{Y}^{-1}\mathbf{X})_{ij}.

Value

Returns a scalar (the ACY error).

Warnings

If Y is singular, function will produce an error.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Amari, S., Cichocki, A., & Yang, H.H. (1996). A new learning algorithm for blind signal separation. In D. S. Touretzky, M. C. Mozer, and M. E. Hasselmo (Eds.), Advances in Neural Information Processing Systems, 8. Cambridge, MA: MIT Press.

Examples

##########   EXAMPLE   ##########

set.seed(1)
X <- matrix(runif(16),4,4)
Y <- matrix(runif(16),4,4)
Z <- X[,c(3,1,2,4)]%*%diag(1:4)
acy(X,Y)
acy(X,Z)

ICA via FastICA, Infomax, or JADE

Description

Computes ICA decomposition using Hyvarinen's (1999) FastICA algorithm, Bell and Sejnowski's (1995) Information-Maximization (Infomax) algorithm, or Cardoso and Souloumiac's (1993, 1996) Joint Approximate Diagonalization of Eigenmatrices (JADE) algorithm.

Usage

ica(X, nc, method = c("fast", "imax", "jade"), ...)

Arguments

X

Data matrix with n rows (samples) and p columns (variables).

nc

Number of components to extract.

method

Method for decomposition.

...

Additional arguments to be passed to other ICA functions (see Details).

Details

ICA Model The ICA model can be written as X = tcrossprod(S, M) + E, where S contains the source signals, M is the mixing matrix, and E contains the noise signals. Columns of X are assumed to have zero mean. The goal is to find the unmixing matrix W such that columns of S = tcrossprod(X, W) are independent as possible.

Whitening Without loss of generality, we can write M = P %*% R where P is a tall matrix and R is an orthogonal rotation matrix. Letting Q denote the pseudoinverse of P, we can whiten the data using Y = tcrossprod(X, Q). The goal is to find the orthongal rotation matrix R such that the source signal estimates S = Y %*% R are as independent as possible. Note that W = crossprod(R, Q).

Method This is a wrapper function for the functions icafast, icaimax, or icajade. See the corresponding function for details on the method, as well as the available arguments (handled by the ... argument).

Value

S

Matrix of source signal estimates (S = Y %*% R).

M

Estimated mixing matrix.

W

Estimated unmixing matrix (W = crossprod(R, Q)).

Y

Whitened data matrix.

Q

Whitening matrix.

R

Orthogonal rotation matrix.

vafs

Variance-accounted-for by each component.

iter

Number of algorithm iterations.

converged

Logical indicating if algorithm converged.

...

Other arguments (if method = "fast" or method = "imax").

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Bell, A.J. & Sejnowski, T.J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129-1159. doi:10.1162/neco.1995.7.6.1129

Cardoso, J.F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. IEE Proceedings-F, 140(6), 362-370. doi:10.1049/ip-f-2.1993.0054

Cardoso, J.F., & Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization. SIAM Journal on Matrix Analysis and Applications, 17(1), 161-164. doi:10.1137/S0895479893259546

Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009

Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626-634. doi:10.1109/72.761722

See Also

icafast for ICA via FastICA

icaimax for ICA via Infomax

icajade for ICA via JADE

Examples

##########   EXAMPLE 1   ##########

# generate noiseless data (p == r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via different algorithms
imod.fast <- ica(Xmat, nc = 2)
imod.imax <- ica(Xmat, nc = 2, method = "imax")
imod.jade <- ica(Xmat, nc = 2, method = "jade")

# compare mixing matrix recovery
acy(Bmat, imod.fast$M)
acy(Bmat, imod.imax$M)
acy(Bmat, imod.jade$M)

# compare source signal recovery
cor(Amat, imod.fast$S)
cor(Amat, imod.imax$S)
cor(Amat, imod.jade$S)


##########   EXAMPLE 2   ##########

# generate noiseless data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via different algorithms
imod.fast <- ica(Xmat, nc = 2)
imod.imax <- ica(Xmat, nc = 2, method = "imax")
imod.jade <- ica(Xmat, nc = 2, method = "jade")

# compare source signal recovery
cor(Amat, imod.fast$S)
cor(Amat, imod.imax$S)
cor(Amat, imod.jade$S)


##########   EXAMPLE 3   ##########

# generate noisy data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), 100, 2)
Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100)
Xmat <- tcrossprod(Amat,Bmat) + Emat

# ICA via different algorithms
imod.fast <- ica(Xmat, nc = 2)
imod.imax <- ica(Xmat, nc = 2, method = "imax")
imod.jade <- ica(Xmat, nc = 2, method = "jade")

# compare source signal recovery
cor(Amat, imod.fast$S)
cor(Amat, imod.imax$S)
cor(Amat, imod.jade$S)

ICA via FastICA Algorithm

Description

Computes ICA decomposition using Hyvarinen's (1999) FastICA algorithm with various options.

Usage

icafast(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), 
        alg = "par", fun = "logcosh", alpha = 1)

Arguments

X

Data matrix with n rows (samples) and p columns (variables).

nc

Number of components to extract.

center

If TRUE, columns of X are mean-centered before ICA decomposition.

maxit

Maximum number of algorithm iterations to allow.

tol

Convergence tolerance.

Rmat

Initial estimate of the nc-by-nc orthogonal rotation matrix.

alg

Algorithm to use: alg="par" to estimate all nc components in parallel (default) or alg="def" for deflation estimation (i.e., projection pursuit).

fun

Contrast function to use for negentropy approximation: fun="logcosh" for log of hyperbolic cosine, fun="exp" for Gaussian kernel, and fun="kur" for kurtosis.

alpha

Tuning parameter for "logcosh" contrast function (1 <= alpha <= 2).

Details

ICA Model The ICA model can be written as X = tcrossprod(S, M) + E, where S contains the source signals, M is the mixing matrix, and E contains the noise signals. Columns of X are assumed to have zero mean. The goal is to find the unmixing matrix W such that columns of S = tcrossprod(X, W) are independent as possible.

Whitening Without loss of generality, we can write M = P %*% R where P is a tall matrix and R is an orthogonal rotation matrix. Letting Q denote the pseudoinverse of P, we can whiten the data using Y = tcrossprod(X, Q). The goal is to find the orthongal rotation matrix R such that the source signal estimates S = Y %*% R are as independent as possible. Note that W = crossprod(R, Q).

FastICA The FastICA algorithm finds the orthogonal rotation matrix R that (approximately) maximizes the negentropy of the estimated source signals. Negentropy is approximated using

J(s)=[E(G(s))E(G(z))]2J(s) = [ E(G(s)) - E(G(z)) ]^2

where E denotes the expectation, G is the contrast function, and z is a standard normal variable. See Hyvarinen (1999) or Helwig and Hong (2013) for specifics of fixed-point algorithm.

Value

S

Matrix of source signal estimates (S = Y %*% R).

M

Estimated mixing matrix.

W

Estimated unmixing matrix (W = crossprod(R, Q)).

Y

Whitened data matrix.

Q

Whitening matrix.

R

Orthogonal rotation matrix.

vafs

Variance-accounted-for by each component.

iter

Number of algorithm iterations.

alg

Algorithm used (same as input).

fun

Contrast function (same as input).

alpha

Tuning parameter (same as input).

converged

Logical indicating if algorithm converged.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009

Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626-634. doi:10.1109/72.761722

See Also

icaimax for ICA via Infomax

icajade for ICA via JADE

Examples

##########   EXAMPLE 1   ##########

# generate noiseless data (p == r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
acy(Bmat, imod$M)
cor(Amat, imod$S)



##########   EXAMPLE 2   ##########

# generate noiseless data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
cor(Amat, imod$S)



##########   EXAMPLE 3   ##########

# generate noisy data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), 100, 2)
Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100)
Xmat <- tcrossprod(Amat,Bmat) + Emat

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
cor(Amat, imod$S)

ICA via Infomax Algorithm

Description

Computes ICA decomposition using Bell and Sejnowski's (1995) Information-Maximization (Infomax) approach with various options.

Usage

icaimax(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), 
        alg = "newton", fun = "tanh", signs = rep(1, nc), signswitch = TRUE, 
        rate = 1, rateanneal = NULL)

Arguments

X

Data matrix with n rows (samples) and p columns (variables).

nc

Number of components to extract.

center

If TRUE, columns of X are mean-centered before ICA decomposition.

maxit

Maximum number of algorithm iterations to allow.

tol

Convergence tolerance.

Rmat

Initial estimate of the nc-by-nc orthogonal rotation matrix.

alg

Algorithm to use: alg="newton" for Newton iteration, and alg="gradient" for gradient descent.

fun

Nonlinear (squashing) function to use for algorithm: fun="tanh" for hyperbolic tangent, fun="log" for logistic, and fun="ext" for extended Infomax.

signs

Vector of length nc such that signs[j] = 1 if j-th component is super-Gaussian and signs[j] = -1 if j-th component is sub-Gaussian. Only used if fun="ext". Ignored if signswitch=TRUE.

signswitch

If TRUE, the signs vector is automatically determined from the data; otherwise a confirmatory ICA decomposition is calculated using input signs vector. Only used if fun="ext".

rate

Learing rate for gradient descent algorithm. Ignored if alg="newton".

rateanneal

Annealing angle and proportion for gradient descent learing rate (see Details). Ignored if alg="newton".

Details

ICA Model The ICA model can be written as X = tcrossprod(S, M) + E, where S contains the source signals, M is the mixing matrix, and E contains the noise signals. Columns of X are assumed to have zero mean. The goal is to find the unmixing matrix W such that columns of S = tcrossprod(X, W) are independent as possible.

Whitening Without loss of generality, we can write M = P %*% R where P is a tall matrix and R is an orthogonal rotation matrix. Letting Q denote the pseudoinverse of P, we can whiten the data using Y = tcrossprod(X, Q). The goal is to find the orthongal rotation matrix R such that the source signal estimates S = Y %*% R are as independent as possible. Note that W = crossprod(R, Q).

Infomax The Infomax approach finds the orthogonal rotation matrix R that (approximately) maximizes the joint entropy of a nonlinear function of the estimated source signals. See Bell and Sejnowski (1995) and Helwig and Hong (2013) for specifics of algorithms.

Value

S

Matrix of source signal estimates (S = Y %*% R).

M

Estimated mixing matrix.

W

Estimated unmixing matrix (W = crossprod(R, Q)).

Y

Whitened data matrix.

Q

Whitening matrix.

R

Orthogonal rotation matrix.

vafs

Variance-accounted-for by each component.

iter

Number of algorithm iterations.

alg

Algorithm used (same as input).

fun

Contrast function (same as input).

signs

Component signs (same as input).

rate

Learning rate (same as input).

converged

Logical indicating if algorithm converged.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Bell, A.J. & Sejnowski, T.J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129-1159. doi:10.1162/neco.1995.7.6.1129

Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009

See Also

icafast for FastICA

icajade for ICA via JADE

Examples

##########   EXAMPLE 1   ##########

# generate noiseless data (p == r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via Infomax with 2 components
imod <- icaimax(Xmat, nc = 2)
acy(Bmat, imod$M)
cor(Amat, imod$S)



##########   EXAMPLE 2   ##########

# generate noiseless data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via Infomax with 2 components
imod <- icaimax(Xmat, nc = 2)
cor(Amat, imod$S)



##########   EXAMPLE 3   ##########

# generate noisy data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), 100, 2)
Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100)
Xmat <- tcrossprod(Amat,Bmat) + Emat

# ICA via Infomax with 2 components
imod <- icaimax(Xmat, nc = 2)
cor(Amat, imod$S)

ICA via JADE Algorithm

Description

Computes ICA decomposition using Cardoso and Souloumiac's (1993, 1996) Joint Approximate Diagonalization of Eigenmatrices (JADE) approach.

Usage

icajade(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc))

Arguments

X

Data matrix with n rows (samples) and p columns (variables).

nc

Number of components to extract.

center

If TRUE, columns of X are mean-centered before ICA decomposition.

maxit

Maximum number of algorithm iterations to allow.

tol

Convergence tolerance.

Rmat

Initial estimate of the nc-by-nc orthogonal rotation matrix.

Details

ICA Model The ICA model can be written as X = tcrossprod(S, M) + E, where S contains the source signals, M is the mixing matrix, and E contains the noise signals. Columns of X are assumed to have zero mean. The goal is to find the unmixing matrix W such that columns of S = tcrossprod(X, W) are independent as possible.

Whitening Without loss of generality, we can write M = P %*% R where P is a tall matrix and R is an orthogonal rotation matrix. Letting Q denote the pseudoinverse of P, we can whiten the data using Y = tcrossprod(X, Q). The goal is to find the orthongal rotation matrix R such that the source signal estimates S = Y %*% R are as independent as possible. Note that W = crossprod(R, Q).

JADE The JADE approach finds the orthogonal rotation matrix R that (approximately) diagonalizes the cumulant array of the source signals. See Cardoso and Souloumiac (1993,1996) and Helwig and Hong (2013) for specifics of the JADE algorithm.

Value

S

Matrix of source signal estimates (S = Y %*% R).

M

Estimated mixing matrix.

W

Estimated unmixing matrix (W = crossprod(R, Q)).

Y

Whitened data matrix.

Q

Whitening matrix.

R

Orthogonal rotation matrix.

vafs

Variance-accounted-for by each component.

iter

Number of algorithm iterations.

converged

Logical indicating if algorithm converged.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Cardoso, J.F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. IEE Proceedings-F, 140(6), 362-370. doi:10.1049/ip-f-2.1993.0054

Cardoso, J.F., & Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization. SIAM Journal on Matrix Analysis and Applications, 17(1), 161-164. doi:10.1137/S0895479893259546

Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009

See Also

icafast for FastICA

icaimax for ICA via Infomax

Examples

##########   EXAMPLE 1   ##########

# generate noiseless data (p == r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via JADE with 2 components
imod <- icajade(Xmat, nc = 2)
acy(Bmat, imod$M)
cor(Amat, imod$S)



##########   EXAMPLE 2   ##########

# generate noiseless data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via JADE with 2 components
imod <- icajade(Xmat, nc = 2)
cor(Amat, imod$S)



##########   EXAMPLE 3   ##########

# generate noisy data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), 100, 2)
Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100)
Xmat <- tcrossprod(Amat,Bmat) + Emat

# ICA via JADE with 2 components
imod <- icajade(Xmat, nc = 2)
cor(Amat, imod$S)

Plot Densities of Source Signal Distributions

Description

Plot density (pdf) and kurtosis for the 18 source signal distributions used in Bach and Jordan (2002); see icasamp for more information.

Usage

icaplot(xseq = seq(-2,2,length.out=500),
        xlab = "", ylab = "", lty = 1, 
        lwd = 1, col = "black", ...)

Arguments

xseq

Sequence of ordered data values for plotting density.

xlab

X-axis label for plot (default is no label).

ylab

Y-axis label for plot (default is no label).

lty

Line type for each density (scalar or vector of length 18).

lwd

Line width for each density (scalar or vector of length 18).

col

Line color for each density (scalar or vector of length 18).

...

Optional inputs for plot.

Value

Produces a plot with NULL return value.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Bach, F.R. (2002). kernel-ica. MATLAB toolbox (ver 1.2) http://www.di.ens.fr/~fbach/kernel-ica/.

Bach, F.R. & Jordan, M.I. (2002). Kernel independent component analysis. Journal of Machine Learning Research, 3, 1-48.

Examples

## Not run: 
##########   EXAMPLE   ##########

quartz(height=9,width=7)
par(mar=c(3,3,3,3))
icaplot()

## End(Not run)

Sample from Various Source Signal Distributions

Description

Sample observations from the 18 source signal distributions used in Bach and Jordan (2002). Can also return density values and kurtosis for each distribution. Use icaplot to plot distributions.

Usage

icasamp(dname, query = c("rnd","pdf","kur"),
        nsamp = NULL, data = NULL)

Arguments

dname

Distribution name: letter "a" through "r" (see Bach & Jordan, 2002).

query

What to return: query="rnd" for random sample, query="pdf" for density values, and query="kur" for kurtosis.

nsamp

Number of observations to sample. Only used if query="rnd".

data

Data values for density evaluation. Only used if query="pdf".

Details

Inspired by usr_distrib.m from Bach's (2002) kernel-ica MATLAB toolbox.

Value

If query="rnd", returns random sample of size nsamp.

If query="pdf", returns density for input data.

If query="kur", returns kurtosis of distribution.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Bach, F.R. (2002). kernel-ica. MATLAB toolbox (ver 1.2) http://www.di.ens.fr/~fbach/kernel-ica/.

Bach, F.R. & Jordan, M.I. (2002). Kernel independent component analysis. Journal of Machine Learning Research, 3, 1-48.

Examples

##########   EXAMPLE   ##########

# sample 1000 observations from distribution "f"
set.seed(123)
mysamp <- icasamp("f","rnd",nsamp=1000)
xr <- range(mysamp)
hist(mysamp,freq=FALSE,ylim=c(0,.8),breaks=sqrt(1000))

# evaluate density of distribution "f"
xseq <- seq(-5,5,length.out=1000)
mypdf <- icasamp("f","pdf",data=xseq)
lines(xseq,mypdf)

# evaluate kurtosis of distribution "f"
icasamp("f","kur")