| Title: | Scientific Computation Using R |
|---|---|
| Description: | A collection of undergraduate level mathematical routines for quantitative work, covering calculus, distribution functions, random variate generation, linear algebra, differential equations and optimization, sized for one semester. The package is best studied with its source open: every routine is short and written to be read. |
| Authors: | Kyun-Seop Bae [aut, cre] |
| Maintainer: | Kyun-Seop Bae <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.5 |
| Built: | 2026-10-10 07:49:25 UTC |
| Source: | https://github.com/ksbae/mathr |
Incomplete beta function with parameter a and b
betai(a, b, x)betai(a, b, x)
a |
shape constant, distribution parameter |
b |
shape constant, distribution parameter |
x |
to which it integrates |
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.
integrated value
Kyun-Seop Bae <[email protected]>
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.
betai(2, 3, 0.4) pbeta(0.4, 2, 3)betai(2, 3, 0.4) pbeta(0.4, 2, 3)
Transforms a binary IEEE 754 format value stored in an array to a decimal value.
Bin2Dec(b)Bin2Dec(b)
b |
a binary value stored in an array with IEEE 754 format |
IEEE 754 defines how a computer represents a floating point number in binary format. This transforms binary format to decimal value.
A decimal value with attribute "Expression"
Kyun-Seop Bae <[email protected]>
IEEE 754
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.5Bin2Dec(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 the input and output of a function be a scalar, a vector, or a matrix.
ChkFx(func, x)ChkFx(func, x)
func |
function to be checked |
x |
given input |
This is to check the types (scalar, vector, matrix) of a function input and output.
CaseNo |
Case number |
minInput |
minimal input type |
minOutput |
minimal output type |
Given |
type of given value |
Repeat |
repeated evaluation of the function |
Gradient |
type of gradient output |
Kyun-Seop Bae <[email protected]>
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))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 of a symmetric matrix
Chol(M)Chol(M)
M |
a symmetric matrix of real values |
Input should be a symmetric matrix of real values. This is one-to-one transformation.
the lower triangular matrix L such that M = L L'
Kyun-Seop Bae <[email protected]>
M <- matrix(c(4, 2, 2, 3), 2, 2) L <- Chol(M) L L %*% t(L) # back to MM <- matrix(c(4, 2, 2, 3), 2, 2) L <- Chol(M) L L %*% t(L) # back to M
Another simple implementation of combination function like choose
Choose(n, r)Choose(n, r)
n |
total count |
r |
count to choose |
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.
combination of choosing r cases from n cases.
Kyun-Seop Bae <[email protected]>
Choose(10, 5) Choose(100, 50) Choose(1000, 500) Choose(2000, 1000)Choose(10, 5) Choose(100, 50) Choose(1000, 500) Choose(2000, 1000)
Convolute two given sequences and return the real part only.
conv(x, y)conv(x, y)
x |
Input sequence |
y |
Response sequence to unit impulse |
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.
Real part of convolution result.
Kyun-Seop Bae <[email protected]>
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 itx <- 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
Convolute two given sequences.
conv0(x, y)conv0(x, y)
x |
Input sequence |
y |
Response sequence to unit impulse |
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.
convolution result.
Kyun-Seop Bae <[email protected]>
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 itx <- 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
Calculate the density of a beta distribution
Dbeta(x, alph, bet)Dbeta(x, alph, bet)
x |
x value for a beta distribution |
alph |
the shape parameter alpha of a beta distribution |
bet |
the shape parameter beta of a beta distribution |
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.
density for x
Kyun-Seop Bae <[email protected]>
Dbeta(0.3, 2, 3) dbeta(0.3, 2, 3)Dbeta(0.3, 2, 3) dbeta(0.3, 2, 3)
Calculate the mass for k of a binomial distribution
Dbinom(k, n, pe)Dbinom(k, n, pe)
k |
count k of wanted event for a binomial distribution |
n |
total number of trials n of a binomial distribution |
pe |
the probability of wanted event for a binomial distribution |
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.
probability for k
Kyun-Seop Bae <[email protected]>
Dbinom(3, 10, 0.3) dbinom(3, 10, 0.3)Dbinom(3, 10, 0.3) dbinom(3, 10, 0.3)
Deconvolute two given sequences measured at different time points.
dc(t1, c1, t2, c2)dc(t1, c1, t2, c2)
t1 |
time points of sequence 1 |
c1 |
measurement values of sequence 1, to be deconvoluted, Response values |
t2 |
time points of sequence 2 |
c2 |
measurement values of sequence 2, input sequence or unit response sequence |
Two measurement sequences with two different time point sequences. The order matters.
Deconvolution result: input sequence or unit response sequence.
Kyun-Seop Bae <[email protected]>
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)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)
Calculate the density of a chi-square distribution
Dchisq(x2, nu)Dchisq(x2, nu)
x2 |
chi-square value for a chi-square distribution |
nu |
degrees of freedom nu of a chi-square distribution |
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.
density for x2
Kyun-Seop Bae <[email protected]>
Dchisq(3, 5) dchisq(3, 5)Dchisq(3, 5) dchisq(3, 5)
Transforms a decimal value to a binary IEEE 754 format value stored in an array.
Dec2Bin(x, Double=TRUE)Dec2Bin(x, Double=TRUE)
x |
a decimal value |
Double |
set this |
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.
An array containing the binary value.
Kyun-Seop Bae <[email protected]>
IEEE 754
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)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)
Deconvolute two given sequences, and return the real part only.
deconv(z, x)deconv(z, x)
z |
Response sequence to be deconvoluted |
x |
Input or Unit response sequence |
The lengths of the two sequences should be the same. The order matters.
Real part of deconvolution result.
Kyun-Seop Bae <[email protected]>
# 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# 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
Deconvolute two given sequences.
deconv0(z, x)deconv0(z, x)
z |
Response sequence to be deconvoluted |
x |
Input or Unit response sequence |
The lengths of the two sequences should be the same. The order matters.
Deconvolution result.
Kyun-Seop Bae <[email protected]>
# 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# 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
From normalized number to decimal number
DENORM(Norm)DENORM(Norm)
Norm |
Normalized number, a vector of length 2. |
Ref. Plauger PJ The Standard C Library. 1992.
Norm[1] * 2^Norm[2]
Kyun-Seop Bae <[email protected]>
DENORM(c(0.625, 4)) # 0.625 * 2^4 DENORM(NORM(pi))DENORM(c(0.625, 4)) # 0.625 * 2^4 DENORM(NORM(pi))
Get first derivative of scalar-input scalar-valued function using Richardson extrapolation 4th order.
Deriv0(fx, x)Deriv0(fx, x)
fx |
scalar valued function |
x |
x point at which derivative is calculated |
This simplest implementation of 4th order Richardson extrapolation for derivative calculation uses a 4x4 matrix.
derivative value
Kyun-Seop Bae <[email protected]>
fx = function(x) x*x Deriv0(fx, 1) # 2fx = function(x) x*x Deriv0(fx, 1) # 2
Get first derivative of scalar-input scalar-valued function using Richardson extrapolation 4th order.
Deriv1(fx, x)Deriv1(fx, x)
fx |
scalar valued function |
x |
x point at which derivative is calculated |
This implementation of 4th order Richardson extrapolation for derivative calculation, uses a vector of length 4.
derivative value
Kyun-Seop Bae <[email protected]>
fx = function(x) x*x Deriv1(fx, 1) # 2fx = function(x) x*x Deriv1(fx, 1) # 2
Get first derivative of vector-input (multi-variable) scalar-valued function using Richardson extrapolation 4th order.
Deriv2(fx, x)Deriv2(fx, x)
fx |
vector-input (multi-variable) scalar-valued function |
x |
vector x at which gradients are calculated |
This implementation of 4th order Richardson extrapolation for derivative calculation, uses 4 vectors.
gradient values
Kyun-Seop Bae <[email protected]>
fx2 = function(x) sum(x*x) Deriv2(fx2, c(1,1)) # 2, 2fx2 = function(x) sum(x*x) Deriv2(fx2, c(1,1)) # 2, 2
Calculate the density of a F distribution
Df(f, nu1, nu2)Df(f, nu1, nu2)
f |
f value for an F distribution |
nu1 |
numerator degrees of freedom nu1 of an F distribution |
nu2 |
denominator degrees of freedom nu2 of an F distribution |
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.
density for f
Kyun-Seop Bae <[email protected]>
Df(2, 4, 9) df(2, 4, 9)Df(2, 4, 9) df(2, 4, 9)
Calculate the density of a gamma distribution
Dgamma(x, alph, bet=1)Dgamma(x, alph, bet=1)
x |
x value for a gamma distribution |
alph |
the shape parameter alpha of a gamma distribution |
bet |
the rate parameter beta of a gamma distribution |
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.
density for x
Kyun-Seop Bae <[email protected]>
Dgamma(2, 3) # bet defaults to 1 Dgamma(2, 3, 2) # bet is a rate dgamma(2, 3, rate=2)Dgamma(2, 3) # bet defaults to 1 Dgamma(2, 3, 2) # bet is a rate dgamma(2, 3, rate=2)
Calculate the density of a log-normal distribution
Dlnorm(x, mu=0, sig=1)Dlnorm(x, mu=0, sig=1)
x |
x value for a log-normal distribution |
mu |
the mean mu of the logarithm for a log-normal distribution |
sig |
the standard deviation sigma of the logarithm for a log-normal distribution |
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.
density for x
Kyun-Seop Bae <[email protected]>
Dlnorm(1) dlnorm(1)Dlnorm(1) dlnorm(1)
Calculate the density of a normal distribution
Dnorm(x, mu=0, sig=1)Dnorm(x, mu=0, sig=1)
x |
x value for a normal distribution |
mu |
the mean mu of a normal distribution |
sig |
the standard deviation sigma of a normal distribution |
It calculates the density for x value of a normal distribution.
sig must be a single number, not a vector.
density for x
Kyun-Seop Bae <[email protected]>
Dnorm(0) Dnorm(1.96, mu=0, sig=1) dnorm(1.96) # base R, for comparisonDnorm(0) Dnorm(1.96, mu=0, sig=1) dnorm(1.96) # base R, for comparison
Calculate the mass for n of a Poisson distribution
Dpois(n, lam)Dpois(n, lam)
n |
count n for a Poisson distribution |
lam |
mean lambda of a Poisson distribution |
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.
mass (probability) for n
Kyun-Seop Bae <[email protected]>
Dpois(3, 4) dpois(3, 4)Dpois(3, 4) dpois(3, 4)
Calculate the density of a t distribution
Dt(t, nu, mu=0, sig=1)Dt(t, nu, mu=0, sig=1)
t |
t value for a t distribution |
nu |
degrees of freedom |
mu |
the location parameter mu of a t distribution |
sig |
the scale parameter sigma of a t distribution |
It calculates the density for t value of a t distribution.
nu and sig must each be a single number, not a vector.
density for t
Kyun-Seop Bae <[email protected]>
Dt(1.5, 7) dt(1.5, 7)Dt(1.5, 7) dt(1.5, 7)
Ranges of an ellipse
ellipRange(center=c(0, 0), radius=c(2, 1), alpha=0)ellipRange(center=c(0, 0), radius=c(2, 1), alpha=0)
center |
coordinate of center point |
radius |
length of long and short axes |
alpha |
rotation angle of long axis |
This returns ranges of x and y axes.
Column vectors of x and y ranges
Kyun-Seop Bae <[email protected]>
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)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)
Plot ellipse and return coordinates.
ellipse(center=c(0, 0), radius=c(2, 1), alpha=0, npoints=100, add=FALSE, ...)ellipse(center=c(0, 0), radius=c(2, 1), alpha=0, npoints=100, add=FALSE, ...)
center |
coordinate of center point |
radius |
length of long and short axes |
alpha |
rotation angle of long axis |
npoints |
number of points used to plot |
add |
whether to plot on an active device |
... |
arguments to be passed to |
This plots an ellipse and returns coordinates of points used to plot.
add must be a single value, not a vector.
coordinates of points used to plot.
Kyun-Seop Bae <[email protected]>
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)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.
erf(x)erf(x)
x |
a real value |
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.
Integration value
Kyun-Seop Bae <[email protected]>
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.
erf(1) 2*pnorm(sqrt(2)) - 1 # base R, for comparisonerf(1) 2*pnorm(sqrt(2)) - 1 # base R, for comparison
Complementary error function.
erfc(x)erfc(x)
x |
a real value |
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.
Kyun-Seop Bae <[email protected]>
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.
erfc(1) 2*pnorm(-sqrt(2)) erfc(20) # still accurate far into the tailerfc(1) 2*pnorm(-sqrt(2)) erfc(20) # still accurate far into the tail
Retained for backward compatibility. Forwards to erfc.
erfccheb(z)erfccheb(z)
z |
a real value |
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.
complementary error function value, the same as erfc(z)
Kyun-Seop Bae <[email protected]>
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.
erfccheb(1) # same as erfc(1)erfccheb(1) # same as erfc(1)
Eigen values and vectors of a Real Symmetric Matrix by Jacobi Algorithm
evJacobi(A)evJacobi(A)
A |
a symmetric matrix of real values |
Input should be a symmetric matrix of real values.
values |
eigen values |
vectors |
eigen vectors |
Kyun-Seop Bae <[email protected]>
A <- matrix(c(2, 1, 1, 2), 2, 2) r <- evJacobi(A) r$values eigen(A)$valuesA <- matrix(c(2, 1, 1, 2), 2, 2) r <- evJacobi(A) r$values eigen(A)$values
Exponential of a real value.
EXP(x)EXP(x)
x |
a single real number to be exponentiated |
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.
Exponential value of a real value.
Kyun-Seop Bae <[email protected]>
EXP(1) exp(1) EXP(-3) sapply(c(1, 2, 3), EXP) # one value at a timeEXP(1) exp(1) EXP(-3) sapply(c(1, 2, 3), EXP) # one value at a time
This gamma function is a generalization of factorial working with real numbers. This function uses Lanczos' approximation and is fairly short.
GAMMA(z)GAMMA(z)
z |
real number |
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.
Gamma function real value
Kyun-Seop Bae <[email protected]>
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)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)
Log gamma function
gammln(xx)gammln(xx)
xx |
a real positive number to be gammln-ed |
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.
Log gamma
Kyun-Seop Bae <[email protected]>
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.
gammln(5) lgamma(5) # log(4!) = log(24)gammln(5) lgamma(5) # log(4!) = log(24)
Incomplete gamma function with parameter a
gammp(a, x)gammp(a, x)
a |
shape constant, distribution parameter |
x |
to which it integrates |
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.
integrated value
Kyun-Seop Bae <[email protected]>
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.
gammp(2, 3) # P(2, 3) pgamma(3, 2)gammp(2, 3) # P(2, 3) pgamma(3, 2)
Incomplete gamma function, Q(a, x) = 1 - P(a, x)
gammq(a, x)gammq(a, x)
a |
shape constant, distribution parameter |
x |
to which it integrates |
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.
integrated value
Kyun-Seop Bae <[email protected]>
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.
gammq(2, 3) # Q(2, 3) = 1 - P(2, 3) pgamma(3, 2, lower.tail=FALSE)gammq(2, 3) # Q(2, 3) = 1 - P(2, 3) pgamma(3, 2, lower.tail=FALSE)
This is the combination of Grad and Hessian
GenD(func, x)GenD(func, x)
func |
vector-input (multi-variable) scalar-valued function |
x |
x point at which gradients are calculated |
Input is the same to Grad and Hessian functions.
gr |
gradient vector |
Hessian |
Hessian matrix |
Kyun-Seop Bae <[email protected]>
GenD(function(x) sum(x^2), c(1, 2)) # gradient, hessian and more GenD(function(x) sum(x^3), 2)GenD(function(x) sum(x^2), c(1, 2)) # gradient, hessian and more GenD(function(x) sum(x^3), 2)
Get the direction vector for the new line search during the minimization.
GetP(Hv, g)GetP(Hv, g)
Hv |
vector of lower triangular part of hessian |
g |
gradient |
The result is the same as solve(H, g).
Direction vector for a new search
Kyun-Seop Bae <[email protected]>
# 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))# 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))
Integrates with given fx, a, and b
GQuad8(fx, a, b)GQuad8(fx, a, b)
fx |
real valued scalar function |
a |
from |
b |
to |
This gives approximate integration value using Gaussian quadrature. It uses 8 points approximation.
Integration value from a to b with function.
Kyun-Seop Bae <[email protected]>
fx = function(x) exp(-x/100) # exact value: 100 * (1 - exp(-0.24)) = 21.33721 GQuad8(fx, 0, 24)fx = function(x) exp(-x/100) # exact value: 100 * (1 - exp(-0.24)) = 21.33721 GQuad8(fx, 0, 24)
Get gradient using Richardson extrapolation 4th order.
Grad(func, x)Grad(func, x)
func |
vector-input (multi-variable) scalar-valued function |
x |
x point at which gradients are calculated |
This is a simplified implementation of 4th order Richardson extrapolation for gradient calculation in numDeriv package.
gradient value
Kyun-Seop Bae <[email protected]>
fx2 = function(x) sum(x*x) Grad(fx2, c(1, 1)) # 2, 2fx2 = function(x) sum(x*x) Grad(fx2, c(1, 1)) # 2, 2
Get hessian using Richardson extrapolation 4th order.
Hessian(fx, x)Hessian(fx, x)
fx |
vector-input (multi-variable) scalar-valued function |
x |
x point at which hessian is calculated |
This is a simplified implementation of 4th order Richardson extrapolation for hessian calculation. See numDeriv package for more detail.
hessian value
Kyun-Seop Bae <[email protected]>
fx2 = function(x) sum(x*x) Hessian(fx2, c(1, 1)) # matrix(c(2, 0, 0, 2), nrow=2)fx2 = function(x) sum(x*x) Hessian(fx2, c(1, 1)) # matrix(c(2, 0, 0, 2), nrow=2)
Inverse of incomplete beta function with parameter a and b
invbetai(p, a, b)invbetai(p, a, b)
p |
integration value between 0 and 1 |
a |
shape constant, distribution parameter |
b |
shape constant, distribution parameter |
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.
inverse of incomplete beta function
Kyun-Seop Bae <[email protected]>
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.
invbetai(0.5, 2, 3) qbeta(0.5, 2, 3)invbetai(0.5, 2, 3) qbeta(0.5, 2, 3)
Inverse of complementary error function.
inverfc(p)inverfc(p)
p |
a positive real integration value usually between 0 and 2 |
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.
inverse of complementary error function value, p = erfc(x)
Kyun-Seop Bae <[email protected]>
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.
inverfc(0.05) erfc(inverfc(0.05)) # back to 0.05 inverfc(1) # 0inverfc(0.05) erfc(inverfc(0.05)) # back to 0.05 inverfc(1) # 0
Inverse of gammp function
invgammp(p, a)invgammp(p, a)
p |
integration value between 0 and 1 |
a |
shape constant, distribution parameter |
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.
inverse of incomplete gamma function
Kyun-Seop Bae <[email protected]>
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.
invgammp(0.5, 3) qgamma(0.5, 3)invgammp(0.5, 3) qgamma(0.5, 3)
LDL' factorization of a symmetric matrix
LDLT(SymMat)LDLT(SymMat)
SymMat |
a symmetric matrix of real values |
Input should be a symmetric matrix of real values.
L |
lower triangular matrix |
D |
the diagonal of D, as a vector |
Kyun-Seop Bae <[email protected]>
M <- matrix(c(4, 2, 2, 3), 2, 2) r <- LDLT(M) r r$L %*% diag(r$D) %*% t(r$L) # back to MM <- matrix(c(4, 2, 2, 3), 2, 2) r <- LDLT(M) r r$L %*% diag(r$D) %*% t(r$L) # back to M
This log gamma function is for larger real values, beyond the range of GAMMA.
This function uses Lanczos' approximation and is fairly short.
LGAMMA(z)LGAMMA(z)
z |
real number |
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.
natural log value of Gamma function
Kyun-Seop Bae <[email protected]>
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)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)
Natural or common log of a real positive value.
LOG(x, DecFlag=FALSE)LOG(x, DecFlag=FALSE)
x |
a real positive number to be log-ed |
DecFlag |
FALSE for natural log, TRUE for common log |
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.
log value of base e or 10.
Kyun-Seop Bae <[email protected]>
LOG(10) log(10) LOG(10, DecFlag=TRUE) # base 10LOG(10) log(10) LOG(10, DecFlag=TRUE) # base 10
Calculate machine epsilon.
MachEps()MachEps()
none
Machine epsilon + 1 is different from 1. But half of it added to one cannot be differentiated from one.
Machine epsilon
Kyun-Seop Bae <[email protected]>
MachEps()MachEps()
Calculate epsilon using x and Start value.
MachEps2(x = 1, Start = 1)MachEps2(x = 1, Start = 1)
x |
value to be used at the equation of x + eps > x |
Start |
Starting epsilon value. Epsilon is halved by each iteration. |
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.
epsilon
Kyun-Seop Bae <[email protected]>
MachEps2(100, 1) (eps2 = MachEps2(1, 0.8)) 1 + eps2 == 1 # FALSE and less than .Machine$double.epsMachEps2(100, 1) (eps2 = MachEps2(1, 0.8)) 1 + eps2 == 1 # FALSE and less than .Machine$double.eps
Calculate meticulous positive epsilon using bisectional search.
MachEps3(MaxIter = 1000)MachEps3(MaxIter = 1000)
MaxIter |
maximum iteration number |
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.
epsilon with "iteration" count attribute
Kyun-Seop Bae <[email protected]>
MachEps3() # Smallest epsilon in this environment. Smaller than .Machine$double.epsMachEps3() # Smallest epsilon in this environment. Smaller than .Machine$double.eps
Calculate meticulous negative epsilon using bisectional search.
MachEps4(MaxIter = 1000)MachEps4(MaxIter = 1000)
MaxIter |
maximum iteration number |
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.
epsilon with "iteration" count attribute
Kyun-Seop Bae <[email protected]>
MachEps4() # different not just sign, but also magnitude with MachEps3()MachEps4() # different not just sign, but also magnitude with MachEps3()
Multiple linear regression.
mlr(y, x.raw, standardize=0, Plot=FALSE)mlr(y, x.raw, standardize=0, Plot=FALSE)
y |
dependent part, a vector |
x.raw |
independent part, a matrix |
standardize |
standardize method, 0 means no standardize |
Plot |
plot or not |
It shows various outputs.
standardize and Plot must each be a single value, not a
vector.
Variable |
variable names |
Estimate |
estimates |
SE |
standard error |
T |
t value |
p-value |
p-value |
Residual |
residuals |
R-Student |
externally studentized residuals |
hat |
hat value |
Cook's D |
Cook's D |
COV-Ratio |
covariance ratio |
DFFITS |
DFFITS |
DFBETAs |
DFBETAs |
Kyun-Seop Bae <[email protected]>
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 columnset.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
Singular value decomposition by Nash
NashSVD(M)NashSVD(M)
M |
a matrix of real values |
This can be used for a least square problem.
Decomposed result
Kyun-Seop Bae <[email protected]>
M <- matrix(c(1, 2, 3, 4, 5, 6), 3, 2) r <- NashSVD(M) r$d # singular values svd(M)$dM <- matrix(c(1, 2, 3, 4, 5, 6), 3, 2) r <- NashSVD(M) r$d # singular values svd(M)$d
From a decimal number to a normalized number
NORM(x)NORM(x)
x |
a decimal number to be normalized |
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.
a vector c(f, n) with x = f * 2^n
Kyun-Seop Bae <[email protected]>
NORM(10) # 10 = 0.625 * 2^4 DENORM(NORM(10))NORM(10) # 10 = 0.625 * 2^4 DENORM(NORM(10))
Shows how a quasi-Newton type optimization algorithm works.
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)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)
x0 |
initial value |
func |
function to be minimized |
TypF |
typical function value, a scalar, to measure approximate magnitude |
GradTol |
gradient tolerance |
StepTol |
x direction or step tolerance |
FnTol |
y or function tolerance |
ItnLimit |
iteration count limit |
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.
par |
final x vector values |
value |
function value at the final x vector |
FnCount |
function evaluation count |
convergence |
convergence is met or not |
grad |
gradient at the final x vector |
hessian |
hessian at the final x vector |
RelGrad |
relative gradient |
RelStep |
relative step |
DelF |
delta F, function value |
f |
function values |
x |
x vector values |
Kyun-Seop Bae <[email protected]>
Optim0(c(-1.2, 1), function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2)Optim0(c(-1.2, 1), function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2)
Calculate the outer product of two three-dimensional vectors
OuterProd(a, b)OuterProd(a, b)
a |
the first vector |
b |
the second vector |
Outer product is not commutative.
Outer product of two three-dimensional vectors
Kyun-Seop Bae <[email protected]>
# 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))# 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))
Calculate the cumulative probability for x of a beta distribution
Pbeta(x, alph, bet)Pbeta(x, alph, bet)
x |
x value for a beta distribution |
alph |
the shape parameter alpha of a beta distribution |
bet |
the shape parameter beta of a beta distribution |
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.
cumulative probability till x
Kyun-Seop Bae <[email protected]>
Pbeta(0.3, 2, 3) pbeta(0.3, 2, 3)Pbeta(0.3, 2, 3) pbeta(0.3, 2, 3)
Calculate the cumulative probability P(X <= k) of a binomial distribution
Pbinom(k, n, pe)Pbinom(k, n, pe)
k |
count k of wanted event for a binomial distribution |
n |
total number of trials n of a binomial distribution |
pe |
the probability of wanted event for a binomial distribution |
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.
cumulative probability till k
Kyun-Seop Bae <[email protected]>
Pbinom(3, 20, 0.3) # P(X <= 3) pbinom(3, 20, 0.3)Pbinom(3, 20, 0.3) # P(X <= 3) pbinom(3, 20, 0.3)
Calculate the cumulative probability for x2 of a chi-square distribution
Pchisq(x2, nu)Pchisq(x2, nu)
x2 |
chi-square value for a chi-square distribution |
nu |
degrees of freedom nu of a chi-square distribution |
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.
cumulative probability till x2
Kyun-Seop Bae <[email protected]>
Pchisq(3.84, 1) # about 0.95 pchisq(3.84, 1)Pchisq(3.84, 1) # about 0.95 pchisq(3.84, 1)
Calculate the cumulative probability for f of an F distribution
Pf(f, nu1, nu2)Pf(f, nu1, nu2)
f |
f value for an F distribution |
nu1 |
numerator degrees of freedom nu1 of an F distribution |
nu2 |
denominator degrees of freedom nu2 of an F distribution |
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.
cumulative probability till f
Kyun-Seop Bae <[email protected]>
Pf(2, 4, 9) pf(2, 4, 9)Pf(2, 4, 9) pf(2, 4, 9)
Calculate the cumulative probability for x of a gamma distribution
Pgamma(x, alph, bet=1)Pgamma(x, alph, bet=1)
x |
x value for a gamma distribution |
alph |
the shape parameter alpha of a gamma distribution |
bet |
the rate parameter beta of a gamma distribution |
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.
cumulative probability till x
Kyun-Seop Bae <[email protected]>
Pgamma(2, 3) Pgamma(2, 3, 2) # bet is a rate pgamma(2, 3, rate=2)Pgamma(2, 3) Pgamma(2, 3, 2) # bet is a rate pgamma(2, 3, rate=2)
Calculate the cumulative probability for x of a log-normal distribution
Plnorm(x, mu=0, sig=1)Plnorm(x, mu=0, sig=1)
x |
x value for a log-normal distribution |
mu |
the mean mu of the logarithm for a log-normal distribution |
sig |
the standard deviation sigma of the logarithm for a log-normal distribution |
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.
cumulative probability till x
Kyun-Seop Bae <[email protected]>
Plnorm(1) plnorm(1)Plnorm(1) plnorm(1)
Calculate the cumulative probability for x of a normal distribution
Pnorm(x, mu=0, sig=1)Pnorm(x, mu=0, sig=1)
x |
x value for a normal distribution |
mu |
the mean mu of a normal distribution |
sig |
the standard deviation sigma of a normal distribution |
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.
cumulative probability till x
Kyun-Seop Bae <[email protected]>
Pnorm(1.96) pnorm(1.96) Pnorm(100, mu=110, sig=15) # IQ below 100Pnorm(1.96) pnorm(1.96) Pnorm(100, mu=110, sig=15) # IQ below 100
It approximates pnorm using a polynomial equation.
PolyNom3(x)PolyNom3(x)
x |
Z value |
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.
Probability of Z < x
Kyun-Seop Bae <[email protected]>
format(pnorm(2), digits=22) format(PolyNom3(2), digits=22) # Compare with the aboveformat(pnorm(2), digits=22) format(PolyNom3(2), digits=22) # Compare with the above
Calculate the cumulative probability P(X <= n) of a Poisson distribution
Ppois(n, lam)Ppois(n, lam)
n |
count n for a Poisson distribution |
lam |
mean lambda of a Poisson distribution |
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.
cumulative probability till n
Kyun-Seop Bae <[email protected]>
Ppois(0, 4) # P(X <= 0) = exp(-4) Ppois(3, 4) ppois(3, 4)Ppois(0, 4) # P(X <= 0) = exp(-4) Ppois(3, 4) ppois(3, 4)
Calculate the cumulative probability for t of a t distribution
Pt(t, nu, mu=0, sig=1)Pt(t, nu, mu=0, sig=1)
t |
t value for a t distribution |
nu |
degrees of freedom |
mu |
the location parameter mu of a t distribution |
sig |
the scale parameter sigma of a t distribution |
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.
cumulative probability till t
Kyun-Seop Bae <[email protected]>
Pt(1.5, 7) pt(1.5, 7)Pt(1.5, 7) pt(1.5, 7)
Calculate the x value for the cumulative probability p of a beta distribution
Qbeta(p, alph, bet)Qbeta(p, alph, bet)
p |
cumulative probability p for a beta distribution |
alph |
the shape parameter alpha of a beta distribution |
bet |
the shape parameter beta of a beta distribution |
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.
x of Pbeta(x, alph, bet)
Kyun-Seop Bae <[email protected]>
Qbeta(0.5, 2, 3) qbeta(0.5, 2, 3)Qbeta(0.5, 2, 3) qbeta(0.5, 2, 3)
Calculate the k for the cumulative probability p of a binomial distribution
Qbinom(p, n, pe)Qbinom(p, n, pe)
p |
cumulative probability of wanted event for a binomial distribution |
n |
total number of trials n of a binomial distribution |
pe |
the probability of wanted event for a binomial distribution |
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.
the smallest k whose cumulative probability reaches p
Kyun-Seop Bae <[email protected]>
Qbinom(0.5, 20, 0.3) # smallest k with P(X <= k) >= p qbinom(0.5, 20, 0.3)Qbinom(0.5, 20, 0.3) # smallest k with P(X <= k) >= p qbinom(0.5, 20, 0.3)
Calculate the chi-square value for the cumulative probability p of a chi-square distribution
Qchisq(p, nu)Qchisq(p, nu)
p |
cumulative probability p for a chi-square distribution |
nu |
degrees of freedom nu of a chi-square distribution |
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.
x2 of Pchisq(x2, nu)
Kyun-Seop Bae <[email protected]>
Qchisq(0.95, 1) # about 3.84 qchisq(0.95, 1)Qchisq(0.95, 1) # about 3.84 qchisq(0.95, 1)
Calculate the f for the cumulative probability p of an F distribution
Qf(p, nu1, nu2)Qf(p, nu1, nu2)
p |
cumulative probability for f of an F distribution |
nu1 |
numerator degrees of freedom nu1 of an F distribution |
nu2 |
denominator degrees of freedom nu2 of an F distribution |
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.
f of Pf(f, nu1, nu2)
Kyun-Seop Bae <[email protected]>
Qf(0.95, 4, 9) qf(0.95, 4, 9)Qf(0.95, 4, 9) qf(0.95, 4, 9)
Calculate the x value for the cumulative probability p of a gamma distribution
Qgamma(p, alph, bet=1)Qgamma(p, alph, bet=1)
p |
cumulative probability p for a gamma distribution |
alph |
the shape parameter alpha of a gamma distribution |
bet |
the rate parameter beta of a gamma distribution |
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.
x of Pgamma(x, alph, bet)
Kyun-Seop Bae <[email protected]>
Qgamma(0.5, 3) qgamma(0.5, 3)Qgamma(0.5, 3) qgamma(0.5, 3)
Calculate the x value for the cumulative probability p of a log-normal distribution
Qlnorm(p, mu=0, sig=1)Qlnorm(p, mu=0, sig=1)
p |
cumulative probability p for a log-normal distribution |
mu |
the mean mu of the logarithm for a log-normal distribution |
sig |
the standard deviation sigma of the logarithm for a log-normal distribution |
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.
x of Plnorm(x, mu, sig)
Kyun-Seop Bae <[email protected]>
Qlnorm(0.5) qlnorm(0.5)Qlnorm(0.5) qlnorm(0.5)
Calculate the x value for the cumulative probability p of a normal distribution
Qnorm(p, mu=0, sig=1)Qnorm(p, mu=0, sig=1)
p |
cumulative probability p for a normal distribution |
mu |
the mean mu of a normal distribution |
sig |
the standard deviation sigma of a normal distribution |
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.
x of Pnorm(x, mu, sig)
Kyun-Seop Bae <[email protected]>
Qnorm(0.975) qnorm(0.975) Qnorm(0.5, mu=110, sig=15)Qnorm(0.975) qnorm(0.975) Qnorm(0.5, mu=110, sig=15)
Calculate the n for the cumulative probability p of a Poisson distribution
Qpois(p, lam)Qpois(p, lam)
p |
cumulative probability of a Poisson distribution |
lam |
mean lambda of a Poisson distribution |
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.
the smallest n whose cumulative probability reaches p
Kyun-Seop Bae <[email protected]>
Qpois(0.5, 4) # smallest n with P(X <= n) >= p qpois(0.5, 4)Qpois(0.5, 4) # smallest n with P(X <= n) >= p qpois(0.5, 4)
Calculate the t value for the cumulative probability p of a t distribution
Qt(p, nu, mu=0, sig=1)Qt(p, nu, mu=0, sig=1)
p |
cumulative probability p for a t distribution |
nu |
degrees of freedom |
mu |
the location parameter mu of a t distribution |
sig |
the scale parameter sigma of a t distribution |
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.
t of Pt(t, nu, mu, sig)
Kyun-Seop Bae <[email protected]>
Qt(0.975, 7) qt(0.975, 7)Qt(0.975, 7) qt(0.975, 7)
Generate random numbers from a beta distribution
Rbeta(n, alph, bet)Rbeta(n, alph, bet)
n |
count, how many do you want |
alph |
shape parameter alpha of a beta distribution |
bet |
shape parameter beta of a beta distribution |
It generates random numbers from random numbers of other distributions.
alph and bet must each be a single number, not a vector.
a vector of random numbers
Kyun-Seop Bae <[email protected]>
set.seed(1) Rbeta(5, 2, 5) round(mean(Rbeta(10000, 2, 5)), 3) # near 2/(2+5)set.seed(1) Rbeta(5, 2, 5) round(mean(Rbeta(10000, 2, 5)), 3) # near 2/(2+5)
Generate random numbers from an exponential distribution
Rexp(n, alpha)Rexp(n, alpha)
n |
count, how many do you want |
alpha |
rate parameter alpha of an exponential distribution |
It uses inverse-transformation method.
a vector of random numbers
Kyun-Seop Bae <[email protected]>
set.seed(1) Rexp(5, 2) round(mean(Rexp(10000, 2)), 3) # alpha is a rate, near 0.5set.seed(1) Rexp(5, 2) round(mean(Rexp(10000, 2)), 3) # alpha is a rate, near 0.5
Generate random numbers from a gamma distribution
Rgamma(n, alph, bet)Rgamma(n, alph, bet)
n |
count, how many do you want |
alph |
shape parameter alpha of a gamma distribution |
bet |
scale parameter beta of a gamma distribution |
It uses acceptance-rejection method.
alph must be a single number, not a vector.
a vector of random numbers
Kyun-Seop Bae <[email protected]>
Banks J. Handbook of simulation (1998) pp152-153
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 6set.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
Generate random numbers from a gamma distribution by a simple method
Rgamma0(n, alph, bet)Rgamma0(n, alph, bet)
n |
count, how many do you want |
alph |
shape parameter alpha of a gamma distribution, at least 1 |
bet |
scale parameter beta of a gamma distribution |
It uses acceptance-rejection method.
alph must be a single number, not a vector.
a vector of random numbers
Kyun-Seop Bae <[email protected]>
Banks J. Handbook of simulation (1998) pp152-153
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 shapesset.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
Advance State vector by tau time with Deriv
RK4(State, tau, Deriv)RK4(State, tau, Deriv)
State |
state vector |
tau |
time amount to advance |
Deriv |
Derivative functions with the input of |
Deriv should accept State vector as an input.
A new state vector after time advancement by tau
Kyun-Seop Bae <[email protected]>
# 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# 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
Generate random numbers from a multivariate normal distribution
Rmvn(n, Mu, Cov)Rmvn(n, Mu, Cov)
n |
count, how many do you want |
Mu |
mean vector of multivariate normal distribution |
Cov |
variance-covariance matrix of multivariate normal distribution |
It uses Cholesky decomposition.
n must be a single number, not a vector.
a matrix of random vectors
Kyun-Seop Bae <[email protected]>
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.7set.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
Generate random numbers from a normal distribution
Rnorm(n, mu=0, sigma=1)Rnorm(n, mu=0, sigma=1)
n |
count, how many do you want |
mu |
mean of normal distribution |
sigma |
standard deviation of normal distribution |
It is more efficient with even n.
n must be a single number, not a vector.
a vector of random numbers
Kyun-Seop Bae <[email protected]>
set.seed(1) Rnorm(5) set.seed(1) round(mean(Rnorm(10000, mu=1, sigma=2)), 2) # near 1set.seed(1) Rnorm(5) set.seed(1) round(mean(Rnorm(10000, mu=1, sigma=2)), 2) # near 1
Integrates with given fx, a, b, and N
romb(fx, a, b, N=4)romb(fx, a, b, N=4)
fx |
real valued scalar function |
a |
from |
b |
to |
N |
degree; the larger, the more precise |
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.
Integration value from a to b with function.
Kyun-Seop Bae <[email protected]>
fx = function(x) exp(-x/100) # exact value: 100 * (1 - exp(-0.24)) = 21.33721 romb(fx, 0, 24, 10)fx = function(x) exp(-x/100) # exact value: 100 * (1 - exp(-0.24)) = 21.33721 romb(fx, 0, 24, 10)
Rounding function in R is round and its behaviour is different from other common software. This is an ordinary rounding function.
Round(x, digits=0)Round(x, digits=0)
x |
real number to be rounded |
digits |
integer indicating the number of decimal places. 0 means 10^0=1, 2 means 10^-2=0.01, -1 means 10^1 = 10 etc. |
Ordinary rounding function (round half away from zero). For more details, see Wikipedia.
rounded value
Kyun-Seop Bae <[email protected]>
Round(1:10 - 0.5) Round(450, -2) Round(1.235, 2)Round(1:10 - 0.5) Round(450, -2) Round(1.235, 2)
Cumulative probability of the run count for the run test
run.p(m, n, r)run.p(m, n, r)
m |
count of fewer species (minimum value = 0) |
n |
count of more frequent species (minimum value = 1) |
r |
count of run (minimum value = 1) |
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.
probability of run count to be less than or equal to r with m and n. P(Run count <= r | m, n)
Kyun-Seop Bae <[email protected]>
run.p(5, 5, 4) # P(run count <= 4 | m=5, n=5) run.p(5, 5, 10)run.p(5, 5, 4) # P(run count <= 4 | m=5, n=5) run.p(5, 5, 10)
p-value for run test
Run.test(RES)Run.test(RES)
RES |
sequence vector to be tested, usually residual values |
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.
p-value, if it is larger than 0.5, it returns 1 - p.
Kyun-Seop Bae <[email protected]>
set.seed(1) Run.test(rnorm(30)) # random, so a large p-value Run.test(rep(c(1, -1), 15)) # perfectly alternatingset.seed(1) Run.test(rnorm(30)) # random, so a large p-value Run.test(rep(c(1, -1), 15)) # perfectly alternating
Random number generation distributed as Uniform[0, 1) by linear congruential method
runifLC(n = 1, Seed = 123457, a = 16807, m = 2147483647, k = 0, b = 10)runifLC(n = 1, Seed = 123457, a = 16807, m = 2147483647, k = 0, b = 10)
n |
count of random numbers to generate |
Seed |
Seed number. If abs(Seed) is less than 1, Sys.time() is the seed. |
a |
multiplier |
m |
modulus. It should be a prime number. |
k |
increment |
b |
count of initial values to discard |
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.
Uniform[0, 1) random numbers of n
Kyun-Seop Bae <[email protected]>
runifLC(8) runifLC(8, 1) runifLC(8, 2) runifLC(8, 0) # random seedrunifLC(8) runifLC(8, 1) runifLC(8, 2) runifLC(8, 0) # random seed
Fourier transformation using the definition, therefore slow.
sft(y)sft(y)
y |
a sequence of fixed time interval |
This is slow but easy to understand.
Fourier transformation result
Kyun-Seop Bae <[email protected]>
sft(c(1, 2, 3, 4)) # slow Fourier transform fft(c(1, 2, 3, 4)) # base R, for comparisonsft(c(1, 2, 3, 4)) # slow Fourier transform fft(c(1, 2, 3, 4)) # base R, for comparison
Integrates with given fx, a, b, and n
simps13(fx, a, b, n)simps13(fx, a, b, n)
fx |
real valued scalar function |
a |
from |
b |
to |
n |
split count |
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.
Integration value from a to b with function.
Kyun-Seop Bae <[email protected]>
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)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)
Integrates with given fx, a, b, and n
simps38(fx, a, b, n)simps38(fx, a, b, n)
fx |
real valued scalar function |
a |
from |
b |
to |
n |
split count |
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.
Integration value from a to b with function.
Kyun-Seop Bae <[email protected]>
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)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)
Square-rooting of a real positive number.
SQRT(x)SQRT(x)
x |
a real positive number to be square-rooted |
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.
square-rooted value
Kyun-Seop Bae <[email protected]>
SQRT(2) sqrt(2)SQRT(2) sqrt(2)
Determinant of a symmetric matrix
SymDet(SymMat)SymDet(SymMat)
SymMat |
a symmetric matrix of real values |
Input should be a symmetric matrix of real values.
Determinant of the input
Kyun-Seop Bae <[email protected]>
M <- matrix(c(4, 1, 1, 3), 2, 2) SymDet(M) det(M)M <- matrix(c(4, 1, 1, 3), 2, 2) SymDet(M) det(M)
Inverse of a symmetric matrix
SymInv(SymMat)SymInv(SymMat)
SymMat |
a symmetric matrix of real values |
Input should be a symmetric matrix of real values.
Inverse of the input
Kyun-Seop Bae <[email protected]>
M <- matrix(c(4, 1, 1, 3), 2, 2) SymInv(M) solve(M)M <- matrix(c(4, 1, 1, 3), 2, 2) SymInv(M) solve(M)
Solution of a symmetric matrix linear equation
SymSol(H, g)SymSol(H, g)
H |
a symmetric matrix of real values, such as hessian |
g |
a response vector, such as gradient |
Input should be a symmetric matrix of real values and a real vector
Solution of the linear system with a real symmetric coefficient matrix
Kyun-Seop Bae <[email protected]>
M <- matrix(c(4, 1, 1, 3), 2, 2) SymSol(M, c(1, 2)) solve(M, c(1, 2))M <- matrix(c(4, 1, 1, 3), 2, 2) SymSol(M, c(1, 2)) solve(M, c(1, 2))
This function returns factorial value from the stored table.
tableFactorial(n)tableFactorial(n)
n |
non-negative integer to be calculated |
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.
factorial value
Kyun-Seop Bae <[email protected]>
tableFactorial(0) tableFactorial(170) tableFactorial(171)tableFactorial(0) tableFactorial(170) tableFactorial(171)
Integrates with given two sequences: x_i and y_i
trapez0(x, y)trapez0(x, y)
x |
sequence of x |
y |
sequence of y |
This is the simplest implementation of the trapezoidal rule. The lengths of the two sequences should be the same.
Integration value with given sequences.
Kyun-Seop Bae <[email protected]>
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)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)
Integrates with given fx, a, b, and n
trapez1(fx, a, b, n)trapez1(fx, a, b, n)
fx |
real valued scalar function |
a |
from |
b |
to |
n |
split count |
This gives approximate integration value using trapezoidal rule.
n must be a single number, not a vector.
Integration value from a to b with function.
Kyun-Seop Bae <[email protected]>
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)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)
Shows how a quasi-Newton type optimization algorithm works.
VMmin(x0, func, MaxIter=9999, Tol=1e-4)VMmin(x0, func, MaxIter=9999, Tol=1e-4)
x0 |
initial value |
func |
function to be minimized |
MaxIter |
iteration count limit |
Tol |
tolerance |
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.
par |
final x vector values |
value |
function value at the final x vector |
FnCount |
function evaluation count |
GrCount |
gradient evaluation count |
convergence |
convergence is met or not |
grad |
gradient at the final x vector |
Kyun-Seop Bae <[email protected]>
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 dimensionVMmin(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