# R/characteristicfunctions.r In prob: Elementary Probability on Finite Sample Spaces

#### Documented in cfbetacfbinomcfcauchycfchisqcfexpcffcfgammacfgeomcfhypercflnormcflogiscfnbinomcfnormcfpoiscfsignrankcftcfunifcfweibullcfwilcox

```#  Characteristic functions
#  Released under GPL 2 or greater
#  Copyright January 2009, G. Jay Kerns

cfbeta <- function(t, shape1, shape2, ncp = 0){
if (shape1 <=0 || shape2 <=0)
stop("shape1, shape2 must be positive")
if (identical(all.equal(ncp, 0), TRUE)){
#require(fAsianOptions)
kummerM(1i*t, shape1, shape1 + shape2)
} else {
fr <- function(x) cos(t*x)*dbeta(x, shape1, shape2, ncp)
fi <- function(x) sin(t*x)*dbeta(x, shape1, shape2, ncp)
Rp <- integrate(fr, lower = 0, upper = 1)\$value
Ip <- integrate(fi, lower = 0, upper = 1)\$value
return( Rp + 1i*Ip )
}
}

cfbinom <- function(t, size, prob){
if (size <= 0 )
stop("size must be positive")
if (prob < 0 || prob > 1)
stop("prob must be in [0,1]")
(prob*exp(1i*t) + (1 - prob))^size
}

cfcauchy = function(t, location = 0, scale = 1){
if (scale <= 0 )
stop("scale must be positive")
exp(1i*location*t - scale*abs(t))
}

cfchisq <- function(t, df, ncp = 0){
if (df < 0 || ncp < 0  )
stop("df and ncp must be nonnegative")
exp(1i*ncp*t/(1-2i*t))/(1 - 2i*t)^(df/2)
}

cfexp <- function(t, rate = 1){
cfgamma(t, shape = 1, scale = 1/rate)
}

cff <- function(t, df1, df2, ncp, kmax = 10){
if (df1 <= 0 || df2 <= 0  )
stop("df1 and df2 must be positive")
#require(fAsianOptions)
if( identical(all.equal(ncp, 0), TRUE) ){
gamma((df1+df2)/2) / gamma(df2/2) * kummerU(-1i*df2*t/df1, df1/2, 1 - df2/2)
} else {
exp(-ncp/2)*sum((ncp/2)^(0:kmax)/factorial(0:kmax)* kummerM(-1i*df2*t/df1, df1/2 + 0:kmax, -df2/2))
}
}

cfgamma <- function(t, shape, rate = 1, scale = 1/rate){
if (rate <= 0  || scale <= 0)
stop("rate must be positive")
(1-scale*1i*t)^(-shape)
}

cfgeom <- function(t, prob){
cfnbinom(t, size = 1, prob = prob)
}

cfhyper <- function(t, m, n, k){
if (m < 0 || n < 0 || k < 0)
stop("m, n, k must be positive")
hypergeo::hypergeo(-k, -m, n - k + 1, exp(1i*t))/hypergeo::hypergeo(-k, -m, n - k + 1, 1)
}

cflnorm <- function(t, meanlog = 0, sdlog = 1){
if (sdlog <= 0)
stop("sdlog must be positive")
if (identical(all.equal(t, 0), TRUE)){
return(1+0i)
} else {
t <- t * exp(meanlog)
Rp1 <- integrate(function(y) exp(-log(y/t)^2/2/sdlog^2) * cos(y)/y, lower = 0, upper = t )\$value
Rp2 <- integrate(function(y) exp(-log(y*t)^2/2/sdlog^2) * cos(1/y)/y, lower = 0, upper = 1/t )\$value
Ip1 <- integrate(function(y) exp(-log(y/t)^2/2/sdlog^2) * sin(y)/y, lower = 0, upper = t )\$value
Ip2 <- integrate(function(y) exp(-log(y*t)^2/2/sdlog^2) * sin(1/y)/y, lower = 0, upper = 1/t )\$value
return((Rp1 + Rp2 + 1i*(Ip1 + Ip2))/(sqrt(2*pi) * sdlog))
}
}

cflogis <- function(t, location = 0, scale = 1){
if (scale <= 0)
stop("scale must be positive")
ifelse( identical(all.equal(t, 0), TRUE),
return(1),
return(exp(1i*location)*pi*scale*t/sinh(pi*scale*t)))
}

cfnbinom <- function(t, size, prob, mu){
if (size <= 0 )
stop("size must be positive")
if (prob <= 0 || prob > 1)
stop("prob must be in (0,1]")
if (!missing(mu)) {
if (!missing(prob))
stop("'prob' and 'mu' both specified")
prob <- size/(size+mu)
}
(prob/(1-(1-prob)*exp(1i*t)))^size
}

cfnorm <- function(t, mean = 0, sd = 1){
if (sd <= 0)
stop("sd must be positive")
exp(1i*mean - (sd*t)^2/2)
}

cfpois <- function(t, lambda){
if (lambda <= 0)
stop("lambda must be positive")
exp(lambda*(exp(1i*t) - 1))
}

cfsignrank <- function(t, n){
sum(exp(1i*t*0:((n+1)*n/2)) * dsignrank(0:((n+1)*n/2), n))
}

cft <- function(t, df, ncp){
if(missing(ncp)) ncp <- 0
if (df <= 0)
stop("df must be positive")
if (identical(all.equal(ncp, 0), TRUE)){
ifelse(identical(all.equal(t, 0), TRUE), 1+0i,
as.complex(besselK(sqrt(df)*abs(t), df/2)*(sqrt(df)*abs(t))^(df/2)/( gamma(df/2) * 2^(df/2 - 1) ))
)
} else {
fr <- function(x) cos(t*x)*dt(x, df, ncp)
fi <- function(x) sin(t*x)*dt(x, df, ncp)
Rp <- integrate(fr, lower = -Inf, upper = Inf)\$value
Ip <- integrate(fi, lower = -Inf, upper = Inf)\$value
return(Rp + 1i*Ip)
}
}

cfunif <- function(t, min = 0, max = 1){
if (max < min)
stop("min cannot be greater than max")
ifelse( identical(all.equal(t, 0), TRUE),
1+0i,
(exp(1i*t*max) - exp(1i*t*min))/(1i*t*(max - min)))
}

cfweibull <- function(t, shape, scale = 1){
if (shape <= 0 || scale <= 0)
stop("shape and scale must be positive")
fr <- function(x) cos(t*x)*dweibull(x, shape, scale)
fi <- function(x) sin(t*x)*dweibull(x, shape, scale)
Rp <- integrate(fr, lower = 0, upper = Inf)\$value
Ip <- integrate(fi, lower = 0, upper = Inf)\$value
return( Rp + 1i*Ip )
}

cfwilcox <- function(t, m, n){
sum(exp(1i*t*0:(m*n)) * dwilcox(0:(m*n), m, n))
}
```

## Try the prob package in your browser

Any scripts or data that you put into this service are public.

prob documentation built on Aug. 28, 2018, 1:04 a.m.