Package {mathr}

Binary to Decimal Format

Description

Transforms a binary IEEE 754 format value stored in an array to a decimal value.

Usage

Bin2Dec(b)

Arguments

Details

IEEE 754 defines how a computer represents a floating point number in binary format. This transforms binary format to decimal value.

Value

A decimal value with attribute "Expression"

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

IEEE 754

Examples

Bin2Dec(rep(0, 32))                  # +0.0
Bin2Dec(c(1, rep(0, 31)))            # -0.0
Bin2Dec(c(0, rep(1, 8), rep(0, 23))) # +Inf
Bin2Dec(c(1, rep(1, 8), rep(0, 23))) # -Inf
Bin2Dec(c(0, rep(1, 8), rep(1, 23))) # NaN
Bin2Dec(c(1, rep(1, 8), rep(1, 23))) # NaN
Bin2Dec(c(1, rep(1, 30), 0))         # NaN

Bin2Dec(c(0, 0,1,1,1,1,1,1,1,1,1,1, 0,0,0,0, rep(0,48))) # +1
Bin2Dec(c(1, 0,1,1,1,1,1,1,1,1,1,1, 0,0,0,0, rep(0,48))) # -1
Bin2Dec(c(0, 0,1,1,1,1,1,1,1,1,1,1, 1,0,0,0, rep(0,48))) # 1.5
Bin2Dec(c(0, 1,0,0,0,0,0,0,0,0,0,0, 0,0,0,0, rep(0,48))) # 2
Bin2Dec(c(0, 1,0,0,0,0,0,0,0,0,0,1, 1,0,1,0, rep(0,48))) # 6.5

Check function input and output types

Description

Check the input and output of a function be a scalar, a vector, or a matrix.

Usage

ChkFx(func, x)

Arguments

Details

This is to check the types (scalar, vector, matrix) of a function input and output.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx0 = function(x) {
  if (length(x) > 1) return(NULL)
  return(x*x)
}
 
fx1 = function(x) x*x  # independent each other (elementwise operation only)
fx2 = function(x) sum(x)
fx3 = function(x) c(sum(x), mean(x)) # c(mean(x), sd(x))

fx4 = function(x) {
  if (is.vector(x)) { Result = sum(x)
  } else if (is.matrix(x)) { Result = rowSums(x)
  } else { Result = NULL }
  return(Result)
}

fx5 = function(x) {
  if (is.vector(x)) { Result = c(sum(x), mean(x))
  } else if (is.matrix(x)) { Result = cbind(rowSums(x), rowMeans(x))
  } else { Result = NULL }
  return(Result)
}

ChkFx(fx1, 3)
ChkFx(fx1, 1:3)
ChkFx(fx2, 3)
ChkFx(fx2, 1:3)
ChkFx(fx3, 3)
ChkFx(fx3, 1:3)
ChkFx(fx4, 1:3)
ChkFx(fx4, matrix(1:9, nrow=3))
ChkFx(fx5, 1:3)
ChkFx(fx5, matrix(1:9, nrow=3))

Cholesky Decomposition

Description

Cholesky decomposition of a symmetric matrix

Usage

Chol(M)

Arguments

Details

Input should be a symmetric matrix of real values. This is one-to-one transformation.

Value

the lower triangular matrix L such that M = L L'

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(4, 2, 2, 3), 2, 2)
L <- Chol(M)
L
L %*% t(L)                   # back to M

Combination function

Description

Another simple implementation of combination function like choose

Usage

Choose(n, r)

Arguments

Details

This is a simple version of the combination function that does not overflow easily. This function works a value at a time: n and r must each be a single number. Use sapply for a vector.

Value

combination of choosing r cases from n cases.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Choose(10, 5)
Choose(100, 50)
Choose(1000, 500)
Choose(2000, 1000)

Denormalize

Description

From normalized number to decimal number

Usage

DENORM(Norm)

Arguments

Details

Ref. Plauger PJ The Standard C Library. 1992.

Value

Norm[1] * 2^Norm[2]

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

DENORM(c(0.625, 4))          # 0.625 * 2^4
DENORM(NORM(pi))

Density of beta distribution

Description

Calculate the density of a beta distribution

Usage

Dbeta(x, alph, bet)

Arguments

Details

It calculates the density for x value of a beta distribution. This function works a value at a time: x, alph and bet must each be a single number. Use sapply for a vector.

Value

density for x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dbeta(0.3, 2, 3)
dbeta(0.3, 2, 3)

Density of binomial distribution

Description

Calculate the mass for k of a binomial distribution

Usage

Dbinom(k, n, pe)

Arguments

Details

It calculates the probability for count k of a binomial distribution. This function works a value at a time: k, n and pe must each be a single number. Use sapply for a vector.

Value

probability for k

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dbinom(3, 10, 0.3)
dbinom(3, 10, 0.3)

Density of chi-square distribution

Description

Calculate the density of a chi-square distribution

Usage

Dchisq(x2, nu)

Arguments

Details

It calculates the density for chi-square value of a chi-square distribution. This function works a value at a time: x2 and nu must each be a single number. Use sapply for a vector.

Value

density for x2

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dchisq(3, 5)
dchisq(3, 5)

Decimal to Binary Format Vector

Description

Transforms a decimal value to a binary IEEE 754 format value stored in an array.

Usage

Dec2Bin(x, Double=TRUE)

Arguments

Details

IEEE 754 defines how a computer represents a floating point number in binary format. This transforms a decimal value to a binary format array. This function works a value at a time: x and Double must each be a single value. Use sapply for a vector.

Value

An array containing the binary value.

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

IEEE 754

Examples

Dec2Bin(1)
Dec2Bin(-1)
Dec2Bin(1.5)
Dec2Bin(2)
Dec2Bin(6.5)

Dec2Bin(1, FALSE)
Dec2Bin(-1, FALSE)
Dec2Bin(1.5, FALSE)
Dec2Bin(2, FALSE)
Dec2Bin(6.5, FALSE)

First derivative, version 0

Description

Get first derivative of scalar-input scalar-valued function using Richardson extrapolation 4th order.

Usage

Deriv0(fx, x)

Arguments

Details

This simplest implementation of 4th order Richardson extrapolation for derivative calculation uses a 4x4 matrix.

Value

derivative value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) x*x
Deriv0(fx, 1) # 2

First derivative, version 1

Description

Get first derivative of scalar-input scalar-valued function using Richardson extrapolation 4th order.

Usage

Deriv1(fx, x)

Arguments

Details

This implementation of 4th order Richardson extrapolation for derivative calculation, uses a vector of length 4.

Value

derivative value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) x*x
Deriv1(fx, 1) # 2

First derivative, version 2

Description

Get first derivative of vector-input (multi-variable) scalar-valued function using Richardson extrapolation 4th order.

Usage

Deriv2(fx, x)

Arguments

Details

This implementation of 4th order Richardson extrapolation for derivative calculation, uses 4 vectors.

Value

gradient values

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx2 = function(x) sum(x*x)
Deriv2(fx2, c(1,1)) # 2, 2

Density of F distribution

Description

Calculate the density of a F distribution

Usage

Df(f, nu1, nu2)

Arguments

Details

It calculates the density for f value of an F distribution. This function works a value at a time: f, nu1 and nu2 must each be a single number. Use sapply for a vector.

Value

density for f

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Df(2, 4, 9)
df(2, 4, 9)

Density of gamma distribution

Description

Calculate the density of a gamma distribution

Usage

Dgamma(x, alph, bet=1)

Arguments

Details

It calculates the density for x value of a gamma distribution. This function works a value at a time: x, alph and bet must each be a single number. Use sapply for a vector.

Value

density for x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dgamma(2, 3)        # bet defaults to 1
Dgamma(2, 3, 2)     # bet is a rate
dgamma(2, 3, rate=2)

Density of log-normal distribution

Description

Calculate the density of a log-normal distribution

Usage

Dlnorm(x, mu=0, sig=1)

Arguments

Details

It calculates the density for x value of a log-normal distribution. x and sig must each be a single number, not a vector. Use sapply over x for several values.

Value

density for x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dlnorm(1)
dlnorm(1)

Density of normal distribution

Description

Calculate the density of a normal distribution

Usage

Dnorm(x, mu=0, sig=1)

Arguments

Details

It calculates the density for x value of a normal distribution. sig must be a single number, not a vector.

Value

density for x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dnorm(0)
Dnorm(1.96, mu=0, sig=1)
dnorm(1.96)          # base R, for comparison

Density of Poisson distribution

Description

Calculate the mass for n of a Poisson distribution

Usage

Dpois(n, lam)

Arguments

Details

It calculates the mass for n of a Poisson distribution. This function works a value at a time: n and lam must each be a single number. Use sapply for a vector.

Value

mass (probability) for n

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dpois(3, 4)
dpois(3, 4)

Density of t distribution

Description

Calculate the density of a t distribution

Usage

Dt(t, nu, mu=0, sig=1)

Arguments

Details

It calculates the density for t value of a t distribution. nu and sig must each be a single number, not a vector.

Value

density for t

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Dt(1.5, 7)
dt(1.5, 7)

Exponential

Description

Exponential of a real value.

Usage

EXP(x)

Arguments

Details

Current CPUs support this function at the instruction level, so this implementation is for study rather than for speed. Ref. Plauger PJ The Standard C Library. 1992. This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

Exponential value of a real value.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

EXP(1)
exp(1)
EXP(-3)
sapply(c(1, 2, 3), EXP)   # one value at a time

Gamma Function

Description

This gamma function is a generalization of factorial working with real numbers. This function uses Lanczos' approximation and is fairly short.

Usage

GAMMA(z)

Arguments

Details

If z is a positive integer, GAMMA(z) = (z - 1)!. If z is zero or a negative integer, it returns +Inf. Ref: Lanczos g=7 coefficients as used by the GNU Scientific Library, gsl_sf_lngamma. This function works a value at a time: z must be a single number. Use sapply for a vector.

Value

Gamma function real value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

GAMMA(0)
GAMMA(1)   # 0! = 1
GAMMA(2)   # 1! = 1
GAMMA(171) # 170! = 7.257416e+306
format(GAMMA(0.5), digits=22) # = sqrt(pi) = 1.772453850905516
format(gamma(0.5), digits=22)
format(sqrt(pi), digits=22)
format(GAMMA(1.5), digits=22)
format(gamma(1.5), digits=22)

Integration using Gaussian quadrature

Description

Integrates with given fx, a, and b

Usage

GQuad8(fx, a, b)

Arguments

Details

This gives approximate integration value using Gaussian quadrature. It uses 8 points approximation.

Value

Integration value from a to b with function.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
# exact value: 100 * (1 - exp(-0.24)) = 21.33721
GQuad8(fx, 0, 24)

Get Hessian and Gradient

Description

This is the combination of Grad and Hessian

Usage

GenD(func, x)

Arguments

Details

Input is the same to Grad and Hessian functions.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

GenD(function(x) sum(x^2), c(1, 2))     # gradient, hessian and more
GenD(function(x) sum(x^3), 2)

Get the Direction Vector

Description

Get the direction vector for the new line search during the minimization.

Usage

GetP(Hv, g)

Arguments

Details

The result is the same as solve(H, g).

Value

Direction vector for a new search

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

# Hv is the lower triangle of the hessian, column by column
GetP(c(4, 1, 3), c(1, 2))
solve(matrix(c(4, 1, 1, 3), 2, 2), c(1, 2))

Gradient Calculation

Description

Get gradient using Richardson extrapolation 4th order.

Usage

Grad(func, x)

Arguments

Details

This is a simplified implementation of 4th order Richardson extrapolation for gradient calculation in numDeriv package.

Value

gradient value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx2 = function(x) sum(x*x)
Grad(fx2, c(1, 1)) # 2, 2

Hessian Calculation

Description

Get hessian using Richardson extrapolation 4th order.

Usage

Hessian(fx, x)

Arguments

Details

This is a simplified implementation of 4th order Richardson extrapolation for hessian calculation. See numDeriv package for more detail.

Value

hessian value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx2 = function(x) sum(x*x)
Hessian(fx2, c(1, 1)) # matrix(c(2, 0, 0, 2), nrow=2)

LDL' Factorization

Description

LDL' factorization of a symmetric matrix

Usage

LDLT(SymMat)

Arguments

Details

Input should be a symmetric matrix of real values.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(4, 2, 2, 3), 2, 2)
r <- LDLT(M)
r
r$L %*% diag(r$D) %*% t(r$L)  # back to M

Natural Log Gamma Function

Description

This log gamma function is for larger real values, beyond the range of GAMMA. This function uses Lanczos' approximation and is fairly short.

Usage

LGAMMA(z)

Arguments

Details

If abs(z) is less than or equal to 171, it uses GAMMA function. Ref: Lanczos g=7 coefficients as used by the GNU Scientific Library, gsl_sf_lngamma. This function works a value at a time: z must be a single number. Use sapply for a vector.

Value

natural log value of Gamma function

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

LGAMMA(0)
LGAMMA(1)   # log(0!) = 0
LGAMMA(2)   # log(1!) = 0
LGAMMA(171) # log(170!) = 706.5731
LGAMMA(-171) # Inf for all negative integers, as in base lgamma()
format(LGAMMA(0.5), digits=22) # = log(sqrt(pi)) = 0.5723649429247
format(lgamma(0.5), digits=22)
format(log(sqrt(pi)), digits=22)
format(LGAMMA(1.5), digits=22)
format(lgamma(1.5), digits=22)

Logarithm natural or common

Description

Natural or common log of a real positive value.

Usage

LOG(x, DecFlag=FALSE)

Arguments

Details

Current CPUs support this function at the instruction level, so this implementation is for study rather than for speed. Ref. Plauger PJ The Standard C Library. 1992. This function works a value at a time: x and DecFlag must each be a single value. Use sapply for a vector.

Value

log value of base e or 10.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

LOG(10)
log(10)
LOG(10, DecFlag=TRUE)        # base 10

Get Machine Epsilon

Description

Calculate machine epsilon.

Usage

MachEps()

Arguments

none

Details

Machine epsilon + 1 is different from 1. But half of it added to one cannot be differentiated from one.

Value

Machine epsilon

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

MachEps()

Get Epsilon a little more general than Machine Epsilon

Description

Calculate epsilon using x and Start value.

Usage

MachEps2(x = 1, Start = 1)

Arguments

Details

This finds the smallest epsilon which satisfies x + eps > x by halving epsilon for each iteration. This function works a value at a time: x and Start must each be a single number. Use sapply for a vector.

Value

epsilon

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

MachEps2(100, 1)
(eps2 = MachEps2(1, 0.8))
1 + eps2 == 1 # FALSE and less than .Machine$double.eps

Get Meticulous Positive Epsilon

Description

Calculate meticulous positive epsilon using bisectional search.

Usage

MachEps3(MaxIter = 1000)

Arguments

Details

This finds really the smallest epsilon by bisectional search. This function works a value at a time: MaxIter must be a single number. Use sapply for a vector.

Value

epsilon with "iteration" count attribute

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

MachEps3() # Smallest epsilon in this environment. Smaller than .Machine$double.eps

Get Meticulous Negative Epsilon

Description

Calculate meticulous negative epsilon using bisectional search.

Usage

MachEps4(MaxIter = 1000)

Arguments

Details

This finds really the smallest epsilon by bisectional search. This function works a value at a time: MaxIter must be a single number. Use sapply for a vector.

Value

epsilon with "iteration" count attribute

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

MachEps4() # different not just sign, but also magnitude with MachEps3()

Normalize

Description

From a decimal number to a normalized number

Usage

NORM(x)

Arguments

Details

It returns (f,n) for x = f * 2 ^ n where (0.5 <= f <= 1). Returns a vector of length 2. For a negative x, f is negative and its magnitude is as for -x. x must be finite. Ref. Plauger PJ The Standard C Library. 1992. This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

a vector c(f, n) with x = f * 2^n

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

NORM(10)                     # 10 = 0.625 * 2^4
DENORM(NORM(10))

Singular Value Decomposition by Nash

Description

Singular value decomposition by Nash

Usage

NashSVD(M)

Arguments

Details

This can be used for a least square problem.

Value

Decomposed result

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(1, 2, 3, 4, 5, 6), 3, 2)
r <- NashSVD(M)
r$d                          # singular values
svd(M)$d

Prototype of optimization function

Description

Shows how a quasi-Newton type optimization algorithm works.

Usage

Optim0(x0, func, TypF=1, GradTol=.Machine$double.eps^(1/3), 
              StepTol=.Machine$double.eps^(2/3), 
              FnTol=.Machine$double.eps^(1/3), ItnLimit=100)

Arguments

Details

If the changes of between steps are less than the tolerances, it stops. GradTol, StepTol, FnTol and ItnLimit must each be a single number, not a vector.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Optim0(c(-1.2, 1), function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2)

Outer product of two three-dimensional vectors

Description

Calculate the outer product of two three-dimensional vectors

Usage

OuterProd(a, b)

Arguments

Details

Outer product is not commutative.

Value

Outer product of two three-dimensional vectors

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

# the vector cross product of two vectors in three dimensions
OuterProd(c(1, 0, 0), c(0, 1, 0))
OuterProd(c(1, 2, 3), c(4, 5, 6))

Cumulative probability of beta distribution

Description

Calculate the cumulative probability for x of a beta distribution

Usage

Pbeta(x, alph, bet)

Arguments

Details

It calculates the cumulative probability for x of a beta distribution. This function works a value at a time: x, alph and bet must each be a single number. Use sapply for a vector.

Value

cumulative probability till x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pbeta(0.3, 2, 3)
pbeta(0.3, 2, 3)

Cumulative probability of binomial distribution

Description

Calculate the cumulative probability P(X <= k) of a binomial distribution

Usage

Pbinom(k, n, pe)

Arguments

Details

It calculates the cumulative probability for k of a binomial distribution. This function works a value at a time: k, n and pe must each be a single number. Use sapply for a vector.

Value

cumulative probability till k

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pbinom(3, 20, 0.3)           # P(X <= 3)
pbinom(3, 20, 0.3)

Cumulative probability of chi-square distribution

Description

Calculate the cumulative probability for x2 of a chi-square distribution

Usage

Pchisq(x2, nu)

Arguments

Details

It calculates the cumulative probability for x2 of a chi-square distribution. This function works a value at a time: x2 and nu must each be a single number. Use sapply for a vector.

Value

cumulative probability till x2

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pchisq(3.84, 1)              # about 0.95
pchisq(3.84, 1)

Cumulative probability of F distribution

Description

Calculate the cumulative probability for f of an F distribution

Usage

Pf(f, nu1, nu2)

Arguments

Details

It calculates the cumulative probability for f of an F distribution. This function works a value at a time: f, nu1 and nu2 must each be a single number. Use sapply for a vector.

Value

cumulative probability till f

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pf(2, 4, 9)
pf(2, 4, 9)

Cumulative probability of gamma distribution

Description

Calculate the cumulative probability for x of a gamma distribution

Usage

Pgamma(x, alph, bet=1)

Arguments

Details

It calculates the cumulative probability for x of a gamma distribution. This function works a value at a time: x, alph and bet must each be a single number. Use sapply for a vector.

Value

cumulative probability till x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pgamma(2, 3)
Pgamma(2, 3, 2)              # bet is a rate
pgamma(2, 3, rate=2)

Cumulative probability of log-normal distribution

Description

Calculate the cumulative probability for x of a log-normal distribution

Usage

Plnorm(x, mu=0, sig=1)

Arguments

Details

It calculates the cumulative probability for x of a log-normal distribution. This function works a value at a time: x, mu and sig must each be a single number. Use sapply for a vector.

Value

cumulative probability till x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Plnorm(1)
plnorm(1)

Cumulative probability of normal distribution

Description

Calculate the cumulative probability for x of a normal distribution

Usage

Pnorm(x, mu=0, sig=1)

Arguments

Details

It calculates the cumulative probability for x of a normal distribution. This function works a value at a time: x, mu and sig must each be a single number. Use sapply for a vector.

Value

cumulative probability till x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pnorm(1.96)
pnorm(1.96)
Pnorm(100, mu=110, sig=15)   # IQ below 100

Approximation of pnorm using a polynomial equation

Description

It approximates pnorm using a polynomial equation.

Usage

PolyNom3(x)

Arguments

Details

It calculates the approximation of pnorm(x). This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

Probability of Z < x

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

format(pnorm(2), digits=22)
format(PolyNom3(2), digits=22) # Compare with the above

Cumulative probability of Poisson distribution

Description

Calculate the cumulative probability P(X <= n) of a Poisson distribution

Usage

Ppois(n, lam)

Arguments

Details

It calculates the cumulative probability for n of a Poisson distribution. This function works a value at a time: n and lam must each be a single number. Use sapply for a vector.

Value

cumulative probability till n

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Ppois(0, 4)                  # P(X <= 0) = exp(-4)
Ppois(3, 4)
ppois(3, 4)

Cumulative probability of t distribution

Description

Calculate the cumulative probability for t of a t distribution

Usage

Pt(t, nu, mu=0, sig=1)

Arguments

Details

It calculates the cumulative probability for t of a t distribution. This function works a value at a time: t, nu, mu and sig must each be a single number. Use sapply for a vector.

Value

cumulative probability till t

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Pt(1.5, 7)
pt(1.5, 7)

Reverse of Pbeta function

Description

Calculate the x value for the cumulative probability p of a beta distribution

Usage

Qbeta(p, alph, bet)

Arguments

Details

It calculates the x value for the cumulative probability p of a beta distribution. This function works a value at a time: p, alph and bet must each be a single number. Use sapply for a vector.

Value

x of Pbeta(x, alph, bet)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qbeta(0.5, 2, 3)
qbeta(0.5, 2, 3)

Reverse of Pbinom function

Description

Calculate the k for the cumulative probability p of a binomial distribution

Usage

Qbinom(p, n, pe)

Arguments

Details

It returns the smallest k with Pbinom(k, n, pe) >= p, the quantile convention of stats::qbinom.

Ties are the one place this can part from base R. Feeding a cumulative probability straight back in as p asks which side of that value the answer falls on, and Pbinom is built from the incomplete beta function rather than the saddle point algorithm base R uses, so the two agree only to about 1e-12 there. A tolerance of that size is applied before the comparison. For any p not sitting exactly on a cumulative probability the two agree exactly. This function works a value at a time: p, n and pe must each be a single number. Use sapply for a vector.

Value

the smallest k whose cumulative probability reaches p

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qbinom(0.5, 20, 0.3)         # smallest k with P(X <= k) >= p
qbinom(0.5, 20, 0.3)

Reverse of Pchisq function

Description

Calculate the chi-square value for the cumulative probability p of a chi-square distribution

Usage

Qchisq(p, nu)

Arguments

Details

It calculates the chi-square value for the cumulative probability p of a chi-square distribution. This function works a value at a time: p and nu must each be a single number. Use sapply for a vector.

Value

x2 of Pchisq(x2, nu)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qchisq(0.95, 1)              # about 3.84
qchisq(0.95, 1)

Reverse of Pf function

Description

Calculate the f for the cumulative probability p of an F distribution

Usage

Qf(p, nu1, nu2)

Arguments

Details

It calculates the f for the cumulative probability p of an F distribution. This function works a value at a time: p, nu1 and nu2 must each be a single number. Use sapply for a vector.

Value

f of Pf(f, nu1, nu2)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qf(0.95, 4, 9)
qf(0.95, 4, 9)

Reverse of Pgamma function

Description

Calculate the x value for the cumulative probability p of a gamma distribution

Usage

Qgamma(p, alph, bet=1)

Arguments

Details

It calculates the x value for the cumulative probability p of a gamma distribution. This function works a value at a time: p, alph and bet must each be a single number. Use sapply for a vector.

Value

x of Pgamma(x, alph, bet)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qgamma(0.5, 3)
qgamma(0.5, 3)

Reverse of Plnorm function

Description

Calculate the x value for the cumulative probability p of a log-normal distribution

Usage

Qlnorm(p, mu=0, sig=1)

Arguments

Details

It calculates the x value for the cumulative probability p of a log-normal distribution. p and sig must each be a single number, not a vector. Use sapply over p for several values.

Value

x of Plnorm(x, mu, sig)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qlnorm(0.5)
qlnorm(0.5)

Reverse of Pnorm function

Description

Calculate the x value for the cumulative probability p of a normal distribution

Usage

Qnorm(p, mu=0, sig=1)

Arguments

Details

It calculates the x value for the cumulative probability p of a normal distribution. p and sig must each be a single number, not a vector. Use sapply over p for several values.

Value

x of Pnorm(x, mu, sig)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qnorm(0.975)
qnorm(0.975)
Qnorm(0.5, mu=110, sig=15)

Reverse of Ppois function

Description

Calculate the n for the cumulative probability p of a Poisson distribution

Usage

Qpois(p, lam)

Arguments

Details

It returns the smallest n with Ppois(n, lam) >= p, the quantile convention of stats::qpois.

Ties are the one place this can part from base R. Feeding a cumulative probability straight back in as p asks which side of that value the answer falls on, and Ppois is built from the incomplete gamma function rather than the saddle point algorithm base R uses, so the two agree only to about 1e-12 there. A tolerance of that size is applied before the comparison. For any p not sitting exactly on a cumulative probability the two agree exactly. This function works a value at a time: p and lam must each be a single number. Use sapply for a vector.

Value

the smallest n whose cumulative probability reaches p

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qpois(0.5, 4)                # smallest n with P(X <= n) >= p
qpois(0.5, 4)

Reverse of Pt function

Description

Calculate the t value for the cumulative probability p of a t distribution

Usage

Qt(p, nu, mu=0, sig=1)

Arguments

Details

It calculates the t value for the cumulative probability p of a t distribution. p and nu must each be a single number, not a vector. Use sapply over p for several values.

Value

t of Pt(t, nu, mu, sig)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Qt(0.975, 7)
qt(0.975, 7)

Runge-Kutta 4th Order Routine for Numerical Integration

Description

Advance State vector by tau time with Deriv

Usage

RK4(State, tau, Deriv)

Arguments

Details

Deriv should accept State vector as an input.

Value

A new state vector after time advancement by tau

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

# harmonic oscillator: x'' = -x, state is c(x, x')
Deriv <- function(s) c(s[2], -s[1])
RK4(c(1, 0), 0.1, Deriv)
c(cos(0.1), -sin(0.1))                  # exact

Random deviate from beta distribution

Description

Generate random numbers from a beta distribution

Usage

Rbeta(n, alph, bet)

Arguments

Details

It generates random numbers from random numbers of other distributions. alph and bet must each be a single number, not a vector.

Value

a vector of random numbers

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
Rbeta(5, 2, 5)
round(mean(Rbeta(10000, 2, 5)), 3)            # near 2/(2+5)

Random deviate from exponential distribution

Description

Generate random numbers from an exponential distribution

Usage

Rexp(n, alpha)

Arguments

Details

It uses inverse-transformation method.

Value

a vector of random numbers

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
Rexp(5, 2)
round(mean(Rexp(10000, 2)), 3)                # alpha is a rate, near 0.5

Random deviate from gamma distribution

Description

Generate random numbers from a gamma distribution

Usage

Rgamma(n, alph, bet)

Arguments

Details

It uses acceptance-rejection method. alph must be a single number, not a vector.

Value

a vector of random numbers

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Banks J. Handbook of simulation (1998) pp152-153

Examples

set.seed(1)
Rgamma(5, 3, 2)
# bet is a scale here, so the mean is alph * bet
round(mean(Rgamma(10000, 3, 2)), 2)           # near 6

Random deviate from gamma distribution, simple method

Description

Generate random numbers from a gamma distribution by a simple method

Usage

Rgamma0(n, alph, bet)

Arguments

Details

It uses acceptance-rejection method. alph must be a single number, not a vector.

Value

a vector of random numbers

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Banks J. Handbook of simulation (1998) pp152-153

Examples

set.seed(1)
Rgamma0(5, 3, 2)
round(mean(Rgamma0(10000, 3, 2)), 2)          # near 6
# alph must be at least 1; use Rgamma for smaller shapes

Random deviate from multivariate normal distribution

Description

Generate random numbers from a multivariate normal distribution

Usage

Rmvn(n, Mu, Cov)

Arguments

Details

It uses Cholesky decomposition. n must be a single number, not a vector.

Value

a matrix of random vectors

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
Cov <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
Rmvn(4, c(0, 0), Cov)
round(cor(Rmvn(5000, c(0, 0), Cov))[1, 2], 2)  # near 0.7

Random deviate from normal distribution

Description

Generate random numbers from a normal distribution

Usage

Rnorm(n, mu=0, sigma=1)

Arguments

Details

It is more efficient with even n. n must be a single number, not a vector.

Value

a vector of random numbers

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
Rnorm(5)
set.seed(1)
round(mean(Rnorm(10000, mu=1, sigma=2)), 2)   # near 1

Ordinary Rounding Function

Description

Rounding function in R is round and its behaviour is different from other common software. This is an ordinary rounding function.

Usage

Round(x, digits=0)

Arguments

Details

Ordinary rounding function (round half away from zero). For more details, see Wikipedia.

Value

rounded value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

Round(1:10 - 0.5)
Round(450, -2)
Round(1.235, 2)

Run test on a sequence

Description

p-value for run test

Usage

Run.test(RES)

Arguments

Details

It calculates p-value for run test. Zeros are omitted first. If no nonzero value is left, there are no runs to count and it returns NA.

Value

p-value, if it is larger than 0.5, it returns 1 - p.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
Run.test(rnorm(30))                 # random, so a large p-value
Run.test(rep(c(1, -1), 15))         # perfectly alternating

Square-rooting

Description

Square-rooting of a real positive number.

Usage

SQRT(x)

Arguments

Details

Current CPUs support this function at the instruction level, so this implementation is for study rather than for speed. Ref. Plauger PJ The Standard C Library. 1992. This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

square-rooted value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

SQRT(2)
sqrt(2)

Determinant of a Symmetric Matrix

Description

Determinant of a symmetric matrix

Usage

SymDet(SymMat)

Arguments

Details

Input should be a symmetric matrix of real values.

Value

Determinant of the input

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(4, 1, 1, 3), 2, 2)
SymDet(M)
det(M)

Inverse of a Symmetric Matrix

Description

Inverse of a symmetric matrix

Usage

SymInv(SymMat)

Arguments

Details

Input should be a symmetric matrix of real values.

Value

Inverse of the input

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(4, 1, 1, 3), 2, 2)
SymInv(M)
solve(M)

Solution of a Symmetric Matrix

Description

Solution of a symmetric matrix linear equation

Usage

SymSol(H, g)

Arguments

Details

Input should be a symmetric matrix of real values and a real vector

Value

Solution of the linear system with a real symmetric coefficient matrix

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

M <- matrix(c(4, 1, 1, 3), 2, 2)
SymSol(M, c(1, 2))
solve(M, c(1, 2))

Simplified version of optimization function

Description

Shows how a quasi-Newton type optimization algorithm works.

Usage

VMmin(x0, func, MaxIter=9999, Tol=1e-4)

Arguments

Details

If the changes of between steps are less than the tolerances, it stops. MaxIter and Tol must each be a single number, not a vector.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

VMmin(c(-1.2, 1), function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2)
VMmin(0, function(x) sum((x - 2.5)^2) + 1)   # one dimension

Continued fraction approximation of incomplete beta function

Description

Incomplete beta function with parameter a and b

Usage

betacf(a, b, x)

Arguments

Details

Ref: Abramowitz and Stegun 26.5.8; beta_cont_frac of the GNU Scientific Library, specfunc/beta_inc.c. This function works a value at a time: a, b and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

# I_x(a,b) = x^a (1-x)^b / (a B(a,b)) * betacf(a, b, x)
a <- 2; b <- 3; x <- 0.4
exp(gammln(a+b) - gammln(a) - gammln(b) +
    a*log(x) + b*log(1-x)) * betacf(a, b, x) / a
pbeta(x, a, b)

Incomplete beta function

Description

Incomplete beta function with parameter a and b

Usage

betai(a, b, x)

Arguments

Details

This calls betacf and gammln. Ref: gsl_sf_beta_inc of the GNU Scientific Library, specfunc/beta_inc.c. This function works a value at a time: a, b and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

betai(2, 3, 0.4)
pbeta(0.4, 2, 3)

Incomplete beta function, kept for compatibility

Description

Retained for backward compatibility. Forwards to betai

Usage

betaiapprox(a, b, x)

Arguments

Details

Retained for backward compatibility. The continued fraction of betacf converges over the whole range once the symmetry transform is applied, so large a and b need no method of their own. This function works a value at a time: a, b and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

betaiapprox(2, 3, 0.4)       # same as betai(2, 3, 0.4)

Linear convolution using fft

Description

Convolute two given sequences and return the real part only.

Usage

conv(x, y)

Arguments

Details

One is input and the other is unit response sequence. The lengths of the two sequences should be the same. The order is not important.

Value

Real part of convolution result.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

x <- c(1, 0.5, 0.25, 0, 0, 0)
y <- c(1, 2, 3, 0, 0, 0)
conv(x, y)                          # linear convolution
round(deconv(conv(x, y), x), 10)    # deconv undoes it

Circular convolution using fft

Description

Convolute two given sequences.

Usage

conv0(x, y)

Arguments

Details

One is input and the other is unit response sequence. The lengths of the two sequences should be the same. The order is not important.

Value

convolution result.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

x <- c(1, 0.5, 0.25, 0.125)
y <- c(1, 2, 3, 4)
Re(conv0(x, y))                     # circular convolution, by fft
round(Re(deconv0(conv0(x, y), x)), 10)   # deconv0 undoes it

Deconvolution of two sequences at different time points

Description

Deconvolute two given sequences measured at different time points.

Usage

dc(t1, c1, t2, c2)

Arguments

Details

Two measurement sequences with two different time point sequences. The order matters.

Value

Deconvolution result: input sequence or unit response sequence.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

t1 <- 0:5
c1 <- c(0, 1, 3, 4, 3, 1)           # response
t2 <- 0:5
c2 <- c(0, 1, 1, 0, 0, 0)           # unit response
dc(t1, c1, t2, c2)

Linear deconvolution using fft

Description

Deconvolute two given sequences, and return the real part only.

Usage

deconv(z, x)

Arguments

Details

The lengths of the two sequences should be the same. The order matters.

Value

Real part of deconvolution result.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

# conv truncates the linear convolution back to the input length, so
# pad the tails with zeros to keep the whole of it. Then deconv
# recovers the input exactly.
x <- c(1, 0.5, 0.25, 0, 0, 0)
y <- c(1, 2, 3, 0, 0, 0)
z <- conv(x, y)
z
round(deconv(z, x), 10)             # back to y

Circular deconvolution using fft

Description

Deconvolute two given sequences.

Usage

deconv0(z, x)

Arguments

Details

The lengths of the two sequences should be the same. The order matters.

Value

Deconvolution result.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

# deconv0 is the inverse of conv0: both work on the circular
# convolution, so the round trip is exact.
x <- c(1, 0.5, 0.25, 0.125)
y <- c(1, 2, 3, 4)
z <- conv0(x, y)
round(Re(deconv0(z, x)), 10)        # back to y

Ranges of an Ellipse

Description

Ranges of an ellipse

Usage

ellipRange(center=c(0, 0), radius=c(2, 1), alpha=0)

Arguments

Details

This returns ranges of x and y axes.

Value

Column vectors of x and y ranges

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

ellipRange()
ellipRange(c(0, 0), c(3, 2), pi/4)
ellipRange(c(1, 1), c(3, 2), pi/4)
ellipRange(c(1, 1), c(3, 2), pi/6)

Ellipse

Description

Plot ellipse and return coordinates.

Usage

ellipse(center=c(0, 0), radius=c(2, 1), alpha=0, npoints=100, add=FALSE, ...)

Arguments

Details

This plots an ellipse and returns coordinates of points used to plot. add must be a single value, not a vector.

Value

coordinates of points used to plot.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

ellipse(asp=1)

ellipse(c(0, 0), c(3, 2), pi/4, asp=1)

ellipse(c(1, 1), c(3, 2), pi/4, asp=1)
points(1, 1, pch="+")
abline(v=1, h=1, lty=3)

Error function

Description

Error function.

Usage

erf(x)

Arguments

Details

This calls gammp.

Ref: erf(x) = P(1/2, x^2), the relation the GNU Scientific Library uses between gsl_sf_erf and gsl_sf_gamma_inc_P. This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

Integration value \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

erf(1)
2*pnorm(sqrt(2)) - 1         # base R, for comparison

Complementary error function

Description

Complementary error function.

Usage

erfc(x)

Arguments

Details

1 - erf(x)

This calls gammq.

Ref: erfc(x) = Q(1/2, x^2), the relation the GNU Scientific Library uses between gsl_sf_erfc and gsl_sf_gamma_inc_Q. This function works a value at a time: x must be a single number. Use sapply for a vector.

Value

1 - \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

erfc(1)
2*pnorm(-sqrt(2))
erfc(20)                     # still accurate far into the tail

Complementary error function, kept for compatibility

Description

Retained for backward compatibility. Forwards to erfc.

Usage

erfccheb(z)

Arguments

Details

Retained for backward compatibility. Forwards to erfc, which now uses the incomplete gamma relation over the whole range. This function works a value at a time: z must be a single number. Use sapply for a vector.

Value

complementary error function value, the same as erfc(z)

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

erfccheb(1)                  # same as erfc(1)

Spectral Decomposition of a Real Symmetric Matrix by Jacobi Algorithm

Description

Eigen values and vectors of a Real Symmetric Matrix by Jacobi Algorithm

Usage

evJacobi(A)

Arguments

Details

Input should be a symmetric matrix of real values.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

A <- matrix(c(2, 1, 1, 2), 2, 2)
r <- evJacobi(A)
r$values
eigen(A)$values

Log gamma function

Description

Log gamma function

Usage

gammln(xx)

Arguments

Details

Ref: Lanczos g=7 coefficients as used by the GNU Scientific Library, gsl_sf_lngamma. This function works a value at a time: xx must be a single number. Use sapply for a vector.

Value

Log gamma

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gammln(5)
lgamma(5)                    # log(4!) = log(24)

Regularized lower incomplete gamma function P(a, x)

Description

Incomplete gamma function with parameter a

Usage

gammp(a, x)

Arguments

Details

This calls gser or gcf. Ref: Cephes igam/igamc, as redistributed in ALGLIB specialfunctions.cs. This function works a value at a time: a and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gammp(2, 3)                  # P(2, 3)
pgamma(3, 2)

Incomplete gamma function, kept for compatibility

Description

Retained for backward compatibility. Forwards to gammp or gammq

Usage

gammpapprox(a, x, psig = 1)

Arguments

Details

This calls gammp or gammq. Retained for backward compatibility. The Cephes series and continued fraction behind them converge over the whole range, so a large a needs no method of its own. This function works a value at a time: a, x and psig must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gammpapprox(2, 3, 1)         # same as gammp(2, 3)
gammpapprox(2, 3, 0)         # same as gammq(2, 3)

Regularized upper incomplete gamma function Q(a, x)

Description

Incomplete gamma function, Q(a, x) = 1 - P(a, x)

Usage

gammq(a, x)

Arguments

Details

This calls gser or gcf. Ref: Cephes igam/igamc, as redistributed in ALGLIB specialfunctions.cs. This function works a value at a time: a and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gammq(2, 3)                  # Q(2, 3) = 1 - P(2, 3)
pgamma(3, 2, lower.tail=FALSE)

Continued fraction approximation of incomplete gamma function

Description

Continued fraction approximation of incomplete gamma function, Q(a, x)

Usage

gcf(a, x)

Arguments

Details

This calls gammln. Ref: Abramowitz and Stegun 6.5.31; Cephes igamc, as redistributed in ALGLIB specialfunctions.cs. This function works a value at a time: a and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gcf(2, 5)                    # continued fraction branch, x > a
pgamma(5, 2, lower.tail=FALSE)

Series approximation of incomplete gamma function

Description

Series approximation of incomplete gamma function, P(a, x)

Usage

gser(a, x)

Arguments

Details

This calls gammln. Ref: Abramowitz and Stegun 6.5.29; Cephes igam, as redistributed in ALGLIB specialfunctions.cs. This function works a value at a time: a and x must each be a single number. Use sapply for a vector.

Value

integrated value

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

gser(2, 1)                   # series branch, x <= a
pgamma(1, 2)

Inverse of incomplete beta function

Description

Inverse of incomplete beta function with parameter a and b

Usage

invbetai(p, a, b)

Arguments

Details

This calls gammln and betai. Ref: bracketed Newton steps on betai, following gsl_cdf_beta_Pinv of the GNU Scientific Library. This function works a value at a time: p, a and b must each be a single number. Use sapply for a vector.

Value

inverse of incomplete beta function

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

invbetai(0.5, 2, 3)
qbeta(0.5, 2, 3)

Inverse of complementary error function

Description

Inverse of complementary error function.

Usage

inverfc(p)

Arguments

Details

It goes through the normal quantile, since erfc(x) = 2 Phi(-x sqrt(2)).

Ref: Cephes ndtri, as redistributed in ALGLIB specialfunctions.cs (invnormalcdf). This function works a value at a time: p must be a single number. Use sapply for a vector.

Value

inverse of complementary error function value, p = erfc(x)

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

inverfc(0.05)
erfc(inverfc(0.05))          # back to 0.05
inverfc(1)                   # 0

Inverse of gammp function

Description

Inverse of gammp function

Usage

invgammp(p, a)

Arguments

Details

This calls gammln, gammp, gammq and inverfc. Ref: Wilson and Hilferty (1931) starting value refined by bracketed Newton steps, following gsl_cdf_gamma_Pinv of the GNU Scientific Library. This function works a value at a time: p and a must each be a single number. Use sapply for a vector.

Value

inverse of incomplete gamma function

Author(s)

Kyun-Seop Bae <k@acr.kr>

References

Galassi M, et al. GNU Scientific Library Reference Manual. 3rd ed. Network Theory, 2009. Moshier SL. Methods and Programs for Mathematical Functions. Ellis Horwood, 1989. Abramowitz M, Stegun IA. Handbook of Mathematical Functions. Dover, 1972.

Examples

invgammp(0.5, 3)
qgamma(0.5, 3)

Multiple Linear Regression

Description

Multiple linear regression.

Usage

mlr(y, x.raw, standardize=0, Plot=FALSE)

Arguments

Details

It shows various outputs. standardize and Plot must each be a single value, not a vector.

Value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

set.seed(1)
n <- 30
X <- data.frame(a = rnorm(n, 10, 2), b = rnorm(n, -5, 3))
y <- 1 + 2*X$a - X$b + rnorm(n)
fit <- mlr(y, X)
fit[[1]]                                  # estimates, SE, t, p
coef(summary(lm(y ~ a + b, data = X)))    # base R, for comparison
mlr(y, X, standardize = 1)[[1]]           # centred by column

Integration using Romberg integration

Description

Integrates with given fx, a, b, and N

Usage

romb(fx, a, b, N=4)

Arguments

Details

This gives approximate integration value using Romberg integration. N is the number of rows of the Romberg table; N = 1 is the trapezoidal rule over the whole interval. N must be a single number, not a vector.

Value

Integration value from a to b with function.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
# exact value: 100 * (1 - exp(-0.24)) = 21.33721
romb(fx, 0, 24, 10)

Cumulative probability of run count for run test

Description

Cumulative probability of the run count for the run test

Usage

run.p(m, n, r)

Arguments

Details

It calculates P(Run count <= r | m, n), the cumulative probability of the run count that Run.test turns into a p-value. This function works a value at a time: m, n and r must each be a single number. Use sapply for a vector.

Value

probability of run count to be less than or equal to r with m and n. P(Run count <= r | m, n)

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

run.p(5, 5, 4)                      # P(run count <= 4 | m=5, n=5)
run.p(5, 5, 10)

Random number of uniform distribution [0, 1) by linear congruential method

Description

Random number generation distributed as Uniform[0, 1) by linear congruential method

Usage

  runifLC(n = 1, Seed = 123457, a = 16807, m = 2147483647, k = 0, b = 10)

Arguments

Details

R[i] = X[i]/(m + 1). X[i] = (a*X[i - 1] + k) mod m. Then, period is m - 1. n, Seed, a and b must each be a single number, not a vector. Use sapply over n for several values.

Value

Uniform[0, 1) random numbers of n

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

  runifLC(8)
  runifLC(8, 1)
  runifLC(8, 2)
  runifLC(8, 0) # random seed

Slow Fourier Transformation

Description

Fourier transformation using the definition, therefore slow.

Usage

sft(y)

Arguments

Details

This is slow but easy to understand.

Value

Fourier transformation result

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

sft(c(1, 2, 3, 4))                  # slow Fourier transform
fft(c(1, 2, 3, 4))                  # base R, for comparison

Integration using Simpson's 1/3 formula

Description

Integrates with given fx, a, b, and n

Usage

simps13(fx, a, b, n)

Arguments

Details

This gives approximate integration value using Simpson's 1/3 formula. Argument n should be an even number. n must be a single number, not a vector.

Value

Integration value from a to b with function.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
# exact value: 100 * (1 - exp(-0.24)) = 21.33721
simps13(fx, 0, 24, 10)   # n must be even
simps13(fx, 0, 24, 30)

Integration using Simpson's 3/8 formula

Description

Integrates with given fx, a, b, and n

Usage

simps38(fx, a, b, n)

Arguments

Details

This gives approximate integration value using Simpson's 3/8 formula. Argument n should be a multiple of 3. n must be a single number, not a vector.

Value

Integration value from a to b with function.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
# exact value: 100 * (1 - exp(-0.24)) = 21.33721
simps38(fx, 0, 24, 12)   # n must be a multiple of 3
simps38(fx, 0, 24, 30)

Factorial using Table

Description

This function returns factorial value from the stored table.

Usage

tableFactorial(n)

Arguments

Details

This function returns factorial value from the table known to be precise. This function works a value at a time: n must be a single number. Use sapply for a vector.

Value

factorial value

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

tableFactorial(0)
tableFactorial(170)
tableFactorial(171)

Integration using trapezoidal rule on two sequences

Description

Integrates with given two sequences: x_i and y_i

Usage

trapez0(x, y)

Arguments

Details

This is the simplest implementation of the trapezoidal rule. The lengths of the two sequences should be the same.

Value

Integration value with given sequences.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
x = c(0, 0.5, 1, 2, 3, 4, 6, 8, 12, 24)
y = fx(x)
trapez0(x,y)

Integration using trapezoidal rule on a function

Description

Integrates with given fx, a, b, and n

Usage

trapez1(fx, a, b, n)

Arguments

Details

This gives approximate integration value using trapezoidal rule. n must be a single number, not a vector.

Value

Integration value from a to b with function.

Author(s)

Kyun-Seop Bae <k@acr.kr>

Examples

fx = function(x) exp(-x/100)
# exact value: 100 * (1 - exp(-0.24)) = 21.33721
trapez1(fx, 0, 24, 10)
trapez1(fx, 0, 24, 1000)