| @@ -0,0 +1,71 @@ | ||
| #' Easy conversion of parameters between stabledist and stable | ||
| #' | ||
| #' \code{pm0_to_pm1} has stabledist parameter inputs for pm=0 and returns pm=1 equivalent parameterization. | ||
| #' \code{pm1_to_pm0} has stabledist parameter inputs for pm=1 and returns pm=0 equivalent parameterization. | ||
| #' | ||
| #' | ||
| #' @name Parameter_Conversion_Nolan_pm1_pm0 | ||
| #' @aliases pm1_to_pm0 pm0_to_pm1 | ||
| #' @param a0 the stabledist 'alpha' for pm=0 in 'stabledist' | ||
| #' @param b0 the stabledist 'beta' for pm=0 in 'stabledist' | ||
| #' @param c0 the stabledist 'gamma' for pm=0 in 'stabledist' | ||
| #' @param d0 the stabledist 'delta' for pm=0 in 'stabledist' | ||
| #' @param a1 the stabledist 'alpha' for pm=1 in 'stabledist' | ||
| #' @param b1 the stabledist 'beta' for pm=1 in 'stabledist' | ||
| #' @param c1 the stabledist 'gamma' for pm=1 in 'stabledist' | ||
| #' @param d1 the stabledist 'delta' for pm=1 in 'stabledist' | ||
| #' | ||
| #' @return What you need. See examples. | ||
| #' @export | ||
| #' @examples | ||
| #' q <- -1 | ||
| #' # nolan pm=1 parameters: | ||
| #' a1 <- 1.3 | ||
| #' b1 <- -0.4 | ||
| #' c1 <- 2 | ||
| #' d1 <- 0.75 | ||
| #' # Convert to nolan pm=0 parameters: | ||
| #' pm0 <- pm1_to_pm0(a1,b1,c1,d1) | ||
| #' a0 <- pm0$a0 | ||
| #' b0 <- pm0$b0 | ||
| #' c0 <- pm0$c0 | ||
| #' d0 <- pm0$d0 | ||
| #' # check: | ||
| #' stabledist::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1) | ||
| #' #> [1] 0.1965513 | ||
| #' # only change delta=d0 for pm=0 | ||
| #' stabledist::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0) | ||
| #' stabledist::pstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0) | ||
| #' #> [1] 0.1965513 | ||
| #' stabledist::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1) | ||
| #' #> [1] 0.0572133 | ||
| #' # only change delta=d0 for pm=0 | ||
| #' stabledist::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0) | ||
| #' stabledist::dstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0) | ||
| #' #> [1] 0.0572133 | ||
| pm0_to_pm1 <- function(a0, b0, c0, d0){ | ||
|
|
||
| if(a0!=1){ | ||
| a1 <- a1 | ||
| b1 <- b1 | ||
| c1 <- c1 | ||
| d1 <- d0 - b1*c1*tan(pi*a1/2) | ||
| } | ||
|
|
||
| return(list(a1=a1, b1=b1, c1=c1, d1=d1 )) | ||
|
|
||
| } | ||
| #' @rdname Parameter_Conversion_Nolan_pm1_pm0 | ||
| #' @export pm1_to_pm0 | ||
| pm1_to_pm0 <- function(a1, b1, c1, d1){ | ||
|
|
||
| if(a1!=1){ | ||
| a0 <- a1 | ||
| b0 <- b1 | ||
| c0 <- c1 | ||
| d0 <- d1 + b1*c1*tan(pi*a1/2) | ||
| } | ||
|
|
||
| return(list(a0=a0, b0=b0, c0=c0, d0=d0)) | ||
|
|
||
| } |
| @@ -0,0 +1,100 @@ | ||
| #' Easy conversion of parameters between stabledist and stable | ||
| #' | ||
| #' \code{sd2s} has stabledist parameter inputs and returns stable parameters. | ||
| #' \code{s2sd} has stable parameter inputs and returns stabledist parameters. | ||
| #' | ||
| #' This is a generic function: methods can be defined for it directly | ||
| #' or via the \code{\link{Summary}} group generic. For this to work properly, | ||
| #' the arguments \code{...} should be unnamed, and dispatch is on the | ||
| #' first argument. | ||
| #' | ||
| #' @name Parameter_Conversion | ||
| #' @aliases sd2s s2sd | ||
| #' @param alpha the stabledist 'alpha' | ||
| #' @param beta the stabledist 'beta' | ||
| #' @param gamma the stabledist 'gamma' | ||
| #' @param delta the stabledist 'delta' | ||
| #' @param pm default 1; currently only value supported. the stabledist parameterization 'pm' | ||
| #' @param tail the stable 'tail' analogous to 'alpha' | ||
| #' @param skew the stable 'skew' analogous to 'beta' | ||
| #' @param disp the stable 'disp' analogous to 'gamma' | ||
| #' @param loc the stable 'loc' analogous to 'delta' | ||
| #' | ||
| #' @return What you need. See examples. | ||
| #' @export | ||
| #' @examples | ||
| #' q <- -1 | ||
| #' # nolan pm=1 parameters: | ||
| #' a <- 1.3 | ||
| #' b <- -0.4 | ||
| #' c <- 2 | ||
| #' d <- 0.75 | ||
| #' s <- sd2s(alpha=a, beta=b, gamma=c, delta=d) | ||
| #' stable::pstable(q, tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc) | ||
| #' stabledist::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #' sd <- s2sd(tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc) | ||
| #' stabledist::pstable(q, alpha=sd$alpha, beta=sd$beta , gamma=sd$gamma , delta=sd$delta, pm=1) | ||
| sd2s <- function(alpha, beta, gamma, delta, pm=1){ | ||
|
|
||
| if(pm==1){ | ||
| a=alpha | ||
| b=beta | ||
| c=gamma | ||
| d=delta | ||
|
|
||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- -sign(b)*2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| } | ||
|
|
||
| return(list(tail=a3, skew=b3, disp=c3, loc=d3 )) | ||
|
|
||
| } | ||
| #' @rdname Parameter_Conversion | ||
| #' @export s2sd | ||
| s2sd <- function(tail, skew, disp, loc, pm=1){ | ||
|
|
||
| if(pm==1){ | ||
| #a=alpha | ||
| #b=beta | ||
| #c=gamma | ||
| #d=delta | ||
|
|
||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| #a3 <- a | ||
| #d3 <- d | ||
| ## the others require calcs: | ||
| #DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| #DEL <- sqrt(DEL2) * sign(1-a) | ||
| #eta_a <- min(a, 2-a) | ||
| ## the lindsey-(3) beta: | ||
| #b3 <- -sign(b)*2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| ## the lindsey-(3) scale: | ||
| #c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
|
|
||
| ## reverse! | ||
| a3 <- tail | ||
| b3 <- skew | ||
| c3 <- disp | ||
| d3 <- loc | ||
| ## solve the above for b, c | ||
| eta_a <- min(a3, 2-a3) | ||
| DEL <- cos( pi * a3/2) / cos(b3 * pi * eta_a/2) | ||
| DEL2 <- DEL^2 | ||
| b <- -sign(b3) * sqrt( (DEL2 - cos(pi*a3/2)^2) / sin(pi*a3/2)^2 ) | ||
| c <- c3 * ( cos(pi*a3/2) / DEL )^(1/a3) | ||
| } | ||
|
|
||
| return(list(alpha=a3, beta=b, gamma=c, delta=d3)) | ||
|
|
||
| } |
| @@ -0,0 +1,105 @@ | ||
| #' Stable Distribution | ||
| #' | ||
| #' These functions provide information about the stable distribution with the | ||
| #' location, the dispersion, the skewness and the tail thickness respectively | ||
| #' modelled by the parameters \code{loc}, \code{disp}, \code{skew} and | ||
| #' \code{tail}. | ||
| #' | ||
| #' \code{dstable}, \code{pstable}, \code{qstable} and \code{hstable} compute | ||
| #' the density, the distribution, the quantile and the hazard functions of a | ||
| #' stable variate. \code{rstable} generates random deviates with the prescribed | ||
| #' stable distribution. | ||
| #' | ||
| #' \code{loc} is a location parameter in the same way as the mean in the normal | ||
| #' distribution: it can take any real value. | ||
| #' | ||
| #' \code{disp} is a dispersion parameter in the same way as the standard | ||
| #' deviation in the normal distribution: it can take any positive value. | ||
| #' | ||
| #' \code{skew} is a skewness parameter: it can take any value in \eqn{(-1,1)}. | ||
| #' The distribution is right-skewed, symmetric and left-skewed when \code{skew} | ||
| #' is negative, null or positive respectively. | ||
| #' | ||
| #' \code{tail} is a tail parameter (often named the characteristic exponent): | ||
| #' it can take any value in \eqn{(0,2)} (with \code{tail=1} and \code{tail=2} | ||
| #' yielding the Cauchy and the normal distributions respectively when symmetry | ||
| #' holds). | ||
| #' | ||
| #' If \code{loc}, \code{disp}, \code{skew}, or \code{tail} are not specified | ||
| #' they assume the default values of \eqn{0}, \eqn{1/sqrt(2)}, \eqn{0} and | ||
| #' \eqn{2} respectively. This corresponds to a normal variate with mean\eqn{=0} | ||
| #' and variance\eqn{=1/2 disp^2}. | ||
| #' | ||
| #' The stable characteristic function is given by \deqn{greekphi(t) = i loca t | ||
| #' - disp {|t|}^{tail} [1+i skew sign(t) greekomega(t,tail)]}{phi(t) = i loc t | ||
| #' - disp |t|^tail [1+i skew sign(t) omega(t,tail)]} where | ||
| #' \deqn{greekomega(t,tail) = \frac{2}{\pi} LOG(ABS(t))}{omega(t,tail) = (2/pi) | ||
| #' log|t|} when \code{tail=1}, and \deqn{greekomega(t,tail) = tan(\frac{\pi | ||
| #' tail}{2})}{omega(t,tail) = tan(pi alpha / 2)} otherwise. | ||
| #' | ||
| #' The characteristic function is inverted using Fourier's transform to obtain | ||
| #' the corresponding stable density. This inversion requires the numerical | ||
| #' evaluation of an integral from \eqn{0} to \eqn{\infty}{infinity}. Two | ||
| #' algorithms are proposed for this. The default is Romberg's method | ||
| #' (\code{integration}="Romberg") which is used to evaluate the integral with | ||
| #' an error bounded by \code{eps}. The alternative method is Simpson's | ||
| #' integration (\code{integration}="Simpson"): it approximates the integral | ||
| #' from \eqn{0} to \eqn{\infty}{infinity} by an integral from \eqn{0} to | ||
| #' \code{up} with \code{npt} points subdividing \eqn{(O, up)}. These three | ||
| #' extra arguments -- \code{integration}, \code{up} and \code{npt} -- are only | ||
| #' available when using \code{dstable}. The other functions are all based on | ||
| #' Romberg's algorithm. | ||
| #' | ||
| #' | ||
| #' @aliases dstable pstable qstable hstable rstable | ||
| #' @param y,q vector of quantiles. | ||
| #' @param p vector of probabilites. | ||
| #' @param n number of observations. | ||
| #' @param loc vector of (real) location parameters. | ||
| #' @param disp vector of (positive) dispersion parameters. | ||
| #' @param skew vector of skewness parameters (in [-1,1]). | ||
| #' @param tail vector of parameters (in [0,2]) related to the tail thickness. | ||
| #' @param eps scalar giving the required precision in computation. | ||
| #' @author Philippe Lambert (Catholic University of Louvain, Belgium, | ||
| #' \email{phlambert@@stat.ucl.ac.be}) and Jim Lindsey. | ||
| #' @seealso \code{stablereg} to fit generalized nonlinear | ||
| #' regression models for the stable distribution parameters. | ||
| #' | ||
| #' \code{stable.mode} to compute the mode of a stable | ||
| #' distribution. | ||
| #' @references Lambert, P. and Lindsey, J.K. (1999) Analysing financial returns | ||
| #' using regression models based on non-symmetric stable distributions. Applied | ||
| #' Statistics, 48, 409-424. | ||
| #' @keywords distribution | ||
| #' @examples | ||
| #' | ||
| #' par(mfrow=c(2,2)) | ||
| #' x <- seq(-5,5,by=0.1) | ||
| #' | ||
| #' # Influence of loc (location) | ||
| #' plot(x,dstable(x,loc=-2,disp=1/sqrt(2),skew=-0.8,tail=1.5), | ||
| #' type="l",ylab="",main="Varying LOCation") | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=-0.8,tail=1.5)) | ||
| #' lines(x,dstable(x,loc=2,disp=1/sqrt(2),skew=-0.8,tail=1.5)) | ||
| #' | ||
| #' # Influence of disp (dispersion) | ||
| #' plot(x,dstable(x,loc=0,disp=0.5,skew=0,tail=1.5), | ||
| #' type="l",ylab="",main="Varying DISPersion") | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0,tail=1.5)) | ||
| #' lines(x,dstable(x,loc=0,disp=0.9,skew=0,tail=1.5)) | ||
| #' | ||
| #' # Influence of skew (skewness) | ||
| #' plot(x,dstable(x,loc=0,disp=1/sqrt(2),skew=-0.8,tail=1.5), | ||
| #' type="l",ylab="",main="Varying SKEWness") | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0,tail=1.5)) | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0.8,tail=1.5)) | ||
| #' | ||
| #' # Influence of tail (tail) | ||
| #' plot(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0,tail=0.8), | ||
| #' type="l",ylab="",main="Varying TAIL thickness") | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0,tail=1.5)) | ||
| #' lines(x,dstable(x,loc=0,disp=1/sqrt(2),skew=0,tail=2)) | ||
| #' | ||
|
|
||
|
|
||
|
|
| @@ -0,0 +1,327 @@ | ||
|
|
||
| <!-- README.md is generated from README.Rmd. Please edit README.Rmd --> | ||
| [](https://travis-ci.org/swihart/stable) [](https://cran.r-project.org/package=stable)  | ||
|
|
||
| `stable` R package | ||
| ================== | ||
|
|
||
| This package is intended to be the developmental version to the CRAN version of [Jim Lindsey's stable](http://www.commanster.eu/rcode.html). The .zip files listed on his homepage have been listed as version 1.0 since 2005. For the subsequent maintenance on this github and CRAN, we will start at 1.1. | ||
|
|
||
| To compare this version with the static v1.0 files on [Jim Lindsey's Homepage](http://www.commanster.eu/rcode.html), it may be useful to use [the compare page for this repo's two branches](https://github.com/swihart/stable/compare/jim-lindsey-homepage-version-1.0...master?diff=split&name=master). | ||
|
|
||
| comparisons with `stabledist` R package | ||
| ======================================= | ||
|
|
||
| In brief, the parameters have different names and are transformations for each other. First, the names: | ||
|
|
||
| | stabledist | stable | | ||
| |------------|--------| | ||
| | alpha | tail | | ||
| | beta | skew | | ||
| | gamma | disp | | ||
| | delta | loc | | ||
|
|
||
| If you read the Lindsey PDF in this repo, be aware that location is given the greek letter gamma and scale is given the greek letter delta. The Nolan PDF does the opposite and is used for `stabledist`. | ||
|
|
||
| For some values for some distributions things match up nicely, as we see with Normal and Cauchy: | ||
|
|
||
| normal distribution | ||
| ------------------- | ||
|
|
||
| ``` r | ||
| q <- 3 | ||
| stable::pstable(q, tail =2, skew=0, disp =1, loc =0) | ||
| #> [1] 0.9830526 | ||
| stabledist::pstable(q, alpha=2, beta=0, gamma=1, delta=0) | ||
| #> [1] 0.9830526 | ||
| ``` | ||
|
|
||
| cauchy distribution | ||
| ------------------- | ||
|
|
||
| ``` r | ||
| q <- 3 | ||
| stable::pstable(q, tail =1, skew=0, disp =1, loc =0) | ||
| #> [1] 0.8975836 | ||
| stabledist::pstable(q, alpha=1, beta=0, gamma=1, delta=0) | ||
| #> [1] 0.8975836 | ||
| ``` | ||
|
|
||
| However, to make `stable` equivalent to `stabledist` in general, some transformations are needed. Please see the following examples. Between `stabledist` and `stable`, the `alpha` is equivalent to `tail` and the `delta` is equivalent to `loc` with no transformation. For the `beta` (`skew`) and `gamma` (`disp`) parameters, a transformation is needed to get equivalent calls. Note differences still may exist to numerical accuracy. | ||
|
|
||
| levy cdf | ||
| -------- | ||
|
|
||
| ``` r | ||
| q <- 0.9 | ||
| # nolan pm=1 parameters: | ||
| a <- 0.5 | ||
| b <- 1 | ||
| c <- .25 | ||
| d <- 0.8 | ||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| stable::pstable(q, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 0.1154242 | ||
| stabledist::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> [1] 0.1138462 | ||
| rmutil::plevy(q, m=d, s=c) | ||
| #> [1] 0.1138463 | ||
| # more accuracy!!!!?! | ||
| stable::pstable(q, tail =a, skew=b3, disp =c3, loc =d, eps = 0.13*1e-7) | ||
| #> [1] 0.1138786 | ||
| ``` | ||
|
|
||
| levy pdf | ||
| -------- | ||
|
|
||
| ``` r | ||
| q <- 0.9 | ||
| # nolan pm=1 parameters: | ||
| a <- 0.5 | ||
| b <- 1 | ||
| c <- .25 | ||
| d <- 0.8 | ||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| stable::dstable(q, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 1.806389 | ||
| stabledist::dstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> Warning in uniroot(function(th) log(g(th)), lower = l.th, upper = u.th, : - | ||
| #> Inf replaced by maximally negative value | ||
| #> Warning in uniroot(function(th) log(g(th)), lower = l.th, upper = u.th, : - | ||
| #> Inf replaced by maximally negative value | ||
| #> Warning in .integrate2(g1, lower = a, upper = b, subdivisions = | ||
| #> subdivisions, : roundoff error is detected in the extrapolation table | ||
| #> [1] 1.807224 | ||
| rmutil::dlevy(q, m=d, s=c) | ||
| #> [1] 1.807224 | ||
| ``` | ||
|
|
||
| levy quantile | ||
| ------------- | ||
|
|
||
| ``` r | ||
| p <- .3 | ||
| # nolan pm=1 parameters: | ||
| a <- 0.5 | ||
| b <- 1 | ||
| c <- .25 | ||
| d <- 0.8 | ||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| stable::qstable(p, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 1.031301 | ||
| stabledist::qstable(p, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> [1] 1.032735 | ||
| rmutil::qlevy(p, m=d, s=c) | ||
| #> [1] 1.032733 | ||
| ``` | ||
|
|
||
| play with alpha not 2 and not 1 | ||
| ------------------------------- | ||
|
|
||
| ``` r | ||
| q <- -1.97 | ||
| # nolan pm=1 parameters: | ||
| a <- 0.8 | ||
| b <- 0 | ||
| c <- 1 | ||
| d <- 0 | ||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- 2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| stable::pstable(q, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 0.1722953 | ||
| stabledist::pstable(q, alpha=a, beta=b , gamma=c , delta=d) | ||
| #> [1] 0.1722945 | ||
| ``` | ||
|
|
||
| play with skew | ||
| -------------- | ||
|
|
||
| ``` r | ||
| q <- -1 | ||
| # nolan pm=1 parameters: | ||
| a <- 1.3 | ||
| b <- 0.4 | ||
| c <- 2 | ||
| d <- 0.75 | ||
| # lindsey-(3) page 415 conversion: | ||
| # tail/alpha and location stay the same | ||
| a3 <- a | ||
| d3 <- d | ||
| # the others require calcs: | ||
| DEL2 <- cos(pi/2 * a)^2 + (-b)^2*sin(pi/2 * a)^2 | ||
| DEL <- sqrt(DEL2) * sign(1-a) | ||
| eta_a <- min(a, 2-a) | ||
| # the lindsey-(3) beta: | ||
| b3 <- -sign(b)*2/(pi*eta_a)*acos( cos(pi/2 * a) / DEL ) | ||
| # the lindsey-(3) scale: | ||
| c3 <- ( (DEL*c^a) / cos(pi/2 * a) )^(1/a) | ||
| stable::pstable(q, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 0.4349168 | ||
| stabledist::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> [1] 0.4348957 | ||
| stable::dstable(q, tail =a, skew=b3, disp =c3, loc =d) | ||
| #> [1] 0.1454112 | ||
| stabledist::dstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> [1] 0.1454111 | ||
| ``` | ||
|
|
||
| The example above, but using `sd2s` and `s2sd` | ||
| ---------------------------------------------- | ||
|
|
||
| ``` r | ||
| q <- -1 | ||
| # nolan pm=1 parameters: | ||
| a <- 1.3 | ||
| b <- -0.4 | ||
| c <- 2 | ||
| d <- 0.75 | ||
| # sd2s takes nolan (stabledist) parameters and returns lindsey (stable) | ||
| s <- stable::sd2s(alpha=a, beta=b, gamma=c, delta=d) | ||
| stable::pstable(q, tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc) | ||
| #> [1] 0.196531 | ||
| stabledist::pstable(q, alpha=a, beta=b , gamma=c , delta=d, pm=1) | ||
| #> [1] 0.1965513 | ||
| # s2sd takes lindsey (stable) parameters and returns nolan (stabledist) | ||
| sd <- stable::s2sd(tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc) | ||
| stabledist::pstable(q, alpha=sd$alpha, beta=sd$beta , gamma=sd$gamma , delta=sd$delta, pm=1) | ||
| #> [1] 0.1965513 | ||
| ``` | ||
|
|
||
| pm1\_to\_pm0 | ||
| ------------ | ||
|
|
||
| ``` r | ||
| q <- -1 | ||
| # nolan pm=1 parameters: | ||
| a1 <- 1.3 | ||
| b1 <- -0.4 | ||
| c1 <- 2 | ||
| d1 <- 0.75 | ||
| # for a1 != 1 | ||
| d0 <- d1 + b1*c1*tan(pi*a1/2) | ||
| # Calculate d0 by hand or use pm1_to_pm0(): | ||
| # Convert to nolan pm=0 parameters: | ||
| pm0 <- stable::pm1_to_pm0(a1,b1,c1,d1) | ||
| a0 <- pm0$a0 | ||
| b0 <- pm0$b0 | ||
| c0 <- pm0$c0 | ||
| d0 <- pm0$d0 | ||
| # check: | ||
| stabledist::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1) | ||
| #> [1] 0.1965513 | ||
| # only change delta=d0 for pm=0 | ||
| stabledist::pstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0) | ||
| #> [1] 0.1965513 | ||
| stabledist::pstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0) | ||
| #> [1] 0.1965513 | ||
| stabledist::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d1, pm=1) | ||
| #> [1] 0.0572133 | ||
| # only change delta=d0 for pm=0 | ||
| stabledist::dstable(q, alpha=a1, beta=b1 , gamma=c1 , delta=d0, pm=0) | ||
| #> [1] 0.0572133 | ||
| stabledist::dstable(q, alpha=a0, beta=b0 , gamma=c0 , delta=d0, pm=0) | ||
| #> [1] 0.0572133 | ||
| ``` | ||
|
|
||
| mode of a stable distribution | ||
| ----------------------------- | ||
|
|
||
| ``` r | ||
| q <- -1 | ||
| # nolan pm=1 parameters: | ||
| # a1 <- 1.3 | ||
| # b1 <- 0.4 | ||
| # c1 <- 2 | ||
| # d1 <- 0.75 | ||
| a1 <- 1.3 | ||
| b1 <- .5 | ||
| c1 <- 1 | ||
| d1 <- 0 | ||
| # for a1 != 1 | ||
| d0 <- d1 + b1*c1*tan(pi*a1/2) | ||
| s <- stable::sd2s(alpha=a1, beta=b1, gamma=c1, delta=d1) | ||
| stable::stable.mode(tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc)$ytilde | ||
| #> [1] -1.13224 | ||
| c1*stabledist::stableMode(alpha=a1, beta=b1)+d0 | ||
| #> [1] -1.133257 | ||
| xran <- seq(-2.5,2.6,0.001) | ||
| ysd <- stabledist::dstable(xran, alpha=a1, beta=b1, gamma=c1, delta=d1, pm=1) | ||
| #plot(xran, ysd) | ||
| xran[ysd == max(ysd)] | ||
| #> [1] -1.133 | ||
| ys <- stable::dstable(xran, tail = s$tail, skew=s$skew, disp = s$disp, loc = s$loc) | ||
| #points(xran, ys, col="blue") | ||
| xran[ys == max(ys)] | ||
| #> [1] -1.133 | ||
| ``` |
| @@ -0,0 +1,41 @@ | ||
| # stable R package | ||
| Bruce Swihart | ||
| Feb 2019 | ||
|
|
||
| ## Re-Submission 1 | ||
|
|
||
| * fixed All declared Imports should be used. NOTE. | ||
| * updated stable.r with Roxygen comments | ||
|
|
||
| ## Test environments | ||
| * local OS X install: R version 3.5.2 (2018-12-20) -- "Eggshell Igloo" | ||
| * Ubuntu 14.04.5 LTS (on travis-ci): R version 3.5.2 (2017-01-27) | ||
| * win-builder: R Under development (unstable) (2019-01-31 r76038) | ||
|
|
||
| ## R CMD check results | ||
| There were no ERRORs or WARNINGs or NOTEs. | ||
|
|
||
|
|
||
| ## Downstream dependencies | ||
| There are currently no downstream dependencies for this package. | ||
|
|
||
| # stable R package | ||
| Bruce Swihart | ||
| May 2018 | ||
|
|
||
| ## Submission 1 of v1.1.3: | ||
| Fixed note; added functions. | ||
|
|
||
| ## Test environments | ||
| * local OS X install: R version 3.4.1 (2017-06-30) | ||
| * ubuntu 12.04 (on travis-ci): R version 3.5.0 (2017-01-27) | ||
| * win-builder: R Under development (unstable) (2018-05-20 r74747) -- "Unsuffered Consequences" | ||
|
|
||
| ## R CMD check results | ||
| There were no ERRORs or WARNINGs or Notes. | ||
|
|
||
| ## Downstream dependencies | ||
| There are currently no downstream dependencies for this package. | ||
|
|
||
| ## Miscellaneous | ||
| Added some functions. |
| @@ -0,0 +1 @@ | ||
| q("no") |
| @@ -1,186 +1,202 @@ | ||
| /* | ||
| * stable : A Library of Functions for Stable Distributions | ||
| * Copyright (C) 1998, 1999, 2000, 2001 P. Lambert and J.K. Lindsey | ||
| * | ||
| * This program is free software; you can redistribute it and/or modify | ||
| * it under the terms of the GNU General Public License as published by | ||
| * the Free Software Foundation; either version 2 of the License, or | ||
| * (at your option) any later version. | ||
| * | ||
| * This program is distributed in the hope that it will be useful, | ||
| * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| * GNU General Public License for more details. | ||
| * | ||
| * You should have received a copy of the GNU General Public License | ||
| * along with this program; if not, write to the Free Software | ||
| * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| * | ||
| * SYNOPSIS | ||
| * | ||
| * void stable(int *n, double *y, double *beta, double *alpha, int *npt, | ||
| * double *up, double *eps, int *type, int *err, double *ffy) | ||
| * void pstable(int *n, double *y, double *beta, double *alpha, | ||
| * double *eps, int *err, double *ffy) | ||
| * | ||
| * DESCRIPTION | ||
| * | ||
| * Functions to compute the density and cdf of a stable distribution | ||
| * | ||
| */ | ||
|
|
||
| #include <math.h> | ||
| #include <stddef.h> | ||
| #include "R.h" | ||
|
|
||
| static int nn; | ||
| static double alphai, etai, setai, cetai, yi, yyi; | ||
|
|
||
| static double fcn1(double s){ | ||
| double sa; | ||
| sa=pow(s,alphai); | ||
| return(cos(-yi*s+sa*setai)*exp(-sa*cetai));} | ||
|
|
||
| static double fcn2(double s){ | ||
| double sa; | ||
| sa=pow(s,-alphai); | ||
| return(cos(-yi/s+sa*setai)*exp(-sa*cetai)/(s*s));} | ||
|
|
||
| static double fcn3(double s){ | ||
| double sa; | ||
| sa=pow(s,alphai); | ||
| return((sin(yyi*s-sa*setai)/s)*exp(-sa*cetai));} | ||
|
|
||
| static double fcn4(double s){ | ||
| double sa; | ||
| sa=pow(s,-alphai); | ||
| return((sin(yyi/s-sa*setai)*s)*exp(-sa*cetai)/(s*s));} | ||
|
|
||
| /* integration routines: stripped down version of Romberg integration | ||
| from rmutil library */ | ||
|
|
||
| static void interp(double x[], double fx[], double *f, double *df) | ||
| { | ||
| int i, j, ni=0; | ||
| double diff1, diff2, tmp1, tmp2, lim1, lim2, tab1[5], tab2[5]; | ||
|
|
||
| tmp1=fabs(x[0]); | ||
| for(i=0;i<5;i++){ | ||
| tmp2=fabs(x[i]); | ||
| if(tmp2<tmp1){ | ||
| ni=i; | ||
| tmp1=tmp2;} | ||
| tab1[i]=fx[i]; | ||
| tab2[i]=fx[i];} | ||
| *f=fx[ni--]; | ||
| for(j=0;j<4;j++){ | ||
| for(i=0;i<=4-j;i++){ | ||
| lim1=x[i]; | ||
| lim2=x[i+j+1]; | ||
| diff1=tab1[i+1]-tab2[i]; | ||
| diff2=lim1-lim2; | ||
| if(diff2==0.0)return; | ||
| diff2=diff1/diff2; | ||
| tab2[i]=lim2*diff2; | ||
| tab1[i]=lim1*diff2;} | ||
| *df=2*ni<(2-j)?tab1[ni+1]:tab2[ni--]; | ||
| *f+=*df;}} | ||
|
|
||
| #define FCN(x) ((*fcn)(x)) | ||
|
|
||
| static double evalfn(double (*fcn)(double), int n) | ||
| { | ||
| double x, nn, tmpsum, pnt1, pnt2; | ||
| static double sum; | ||
| int i,j; | ||
|
|
||
| if (n==1){ | ||
| sum=FCN(0.5); | ||
| return(sum);} | ||
| else { | ||
| for(i=1,j=1;j<n-1;j++) i*=3; | ||
| nn=i; | ||
| pnt1=1/(3.0*nn); | ||
| pnt2=2.0*pnt1; | ||
| x=0.5*pnt1; | ||
| tmpsum=0.0; | ||
| for(j=1;j<=i;j++){ | ||
| tmpsum+=FCN(x); | ||
| x+=pnt2; | ||
| tmpsum+=FCN(x); | ||
| x+=pnt1;} | ||
| sum=(sum+tmpsum/nn)/3.0; | ||
| return(sum);}} | ||
|
|
||
| static double romberg(double (*fcn)(double), double eps) | ||
| { | ||
| int j,j1; | ||
| double sum,errsum,x[16],fx[16]; | ||
|
|
||
| x[0]=1.0; | ||
| for(j=0;j<16;j++){ | ||
| j1=j+1; | ||
| fx[j]=evalfn(fcn,j1); | ||
| if(j1>=5){ | ||
| interp(&x[j1-5],&fx[j1-5],&sum,&errsum); | ||
| if(fabs(errsum)<eps*fabs(sum))return(sum);} | ||
| x[j1]=x[j]/9.0; | ||
| fx[j1]=fx[j];} | ||
| sum=R_NaN; | ||
| return(sum);} | ||
|
|
||
| /* density of a stable distribution */ | ||
| void stable(int *n, double *y, double *beta, double *alpha, int *npt, | ||
| double *up, double *eps, int *type, int *err, double *ffy) | ||
| { | ||
| int i, j; | ||
| double h, s, *eta, *seta, *ceta, *sa; | ||
| *err=0; | ||
| eta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| seta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| ceta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| sa=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| nn=*n; | ||
| if(!eta||!seta||!ceta||!sa){ | ||
| *err=1; | ||
| return;} | ||
| for(i=0;i<*n;i++){ | ||
| ffy[i]=0.0; | ||
| eta[i]=beta[i]*(1.0-fabs(1.0-alpha[i]))*M_PI/2.0; | ||
| seta[i]=sin(eta[i]); | ||
| ceta[i]=cos(eta[i]);} | ||
| if(*type==1){ | ||
| *npt=*npt-*npt%2; | ||
| h=*up/ *npt; | ||
| for(j=0;j<=*npt;j++){ | ||
| s=(*npt-j)*h; | ||
| for(i=0;i<*n;i++){ | ||
| sa[i]=pow(s,alpha[i]); | ||
| ffy[i]=ffy[i]+(4-2*(j%2==0)-(j==1||j==*npt))*cos(-y[i]*s+sa[i]*seta[i])*exp(-sa[i]*ceta[i]);}} | ||
| for(i=0;i<*n;i++)ffy[i]=ffy[i]*h/3.0/M_PI;} | ||
| else { | ||
| for(i=0;i<*n;i++){ | ||
| alphai=alpha[i]; | ||
| yi=y[i]; | ||
| setai=seta[i]; | ||
| cetai=ceta[i]; | ||
| ffy[i]=(romberg(fcn1, *eps)+romberg(fcn2, *eps))/M_PI;}}} | ||
|
|
||
| /* cdf of a stable distribution */ | ||
| void pstable(int *n, double *y, double *beta, double *alpha, | ||
| double *eps, int *err, double *ffy) | ||
| { | ||
| int i; | ||
| *err=0; | ||
| nn=*n; | ||
| for(i=0;i<*n;i++){ | ||
| ffy[i]=0.0; | ||
| etai=beta[i]*(1.0-fabs(1.0-alpha[i]))*M_PI/2.0; | ||
| setai=sin(etai); | ||
| cetai=cos(etai); | ||
| alphai=alpha[i]; | ||
| yyi=y[i]; | ||
| if(etai==0.&&yyi==0) | ||
| ffy[i]=0.5; | ||
| else ffy[i]=0.5+(romberg(fcn3, *eps)+romberg(fcn4, *eps))/M_PI;}} | ||
|
|
||
| /* | ||
| * stable : A Library of Functions for Stable Distributions | ||
| * Copyright (C) 1998, 1999, 2000, 2001 P. Lambert and J.K. Lindsey | ||
| * | ||
| * This program is free software; you can redistribute it and/or modify | ||
| * it under the terms of the GNU General Public License as published by | ||
| * the Free Software Foundation; either version 2 of the License, or | ||
| * (at your option) any later version. | ||
| * | ||
| * This program is distributed in the hope that it will be useful, | ||
| * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| * GNU General Public License for more details. | ||
| * | ||
| * You should have received a copy of the GNU General Public License | ||
| * along with this program; if not, write to the Free Software | ||
| * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| * | ||
| * SYNOPSIS | ||
| * | ||
| * void stable(int *n, double *y, double *beta, double *alpha, int *npt, | ||
| * double *up, double *eps, int *type, int *err, double *ffy) | ||
| * void pstable_c(int *n, double *y, double *beta, double *alpha, | ||
| * double *eps, int *err, double *ffy) | ||
| * | ||
| * DESCRIPTION | ||
| * | ||
| * Functions to compute the density and cdf of a stable distribution | ||
| * | ||
| */ | ||
|
|
||
| #include <math.h> | ||
| #include <stddef.h> | ||
| #include "R.h" | ||
|
|
||
| static int nn; | ||
| static double alphai, etai, setai, cetai, yi, yyi; | ||
|
|
||
| static double fcn1(double s){ | ||
| double sa; | ||
| sa=pow(s,alphai); | ||
| return(cos(-yi*s+sa*setai)*exp(-sa*cetai));} | ||
|
|
||
| static double fcn2(double s){ | ||
| double sa; | ||
| sa=pow(s,-alphai); | ||
| return(cos(-yi/s+sa*setai)*exp(-sa*cetai)/(s*s));} | ||
|
|
||
| static double fcn3(double s){ | ||
| double sa; | ||
| sa=pow(s,alphai); | ||
| return((sin(yyi*s-sa*setai)/s)*exp(-sa*cetai));} | ||
|
|
||
| static double fcn4(double s){ | ||
| double sa; | ||
| sa=pow(s,-alphai); | ||
| return((sin(yyi/s-sa*setai)*s)*exp(-sa*cetai)/(s*s));} | ||
|
|
||
| /* integration routines: stripped down version of Romberg integration | ||
| from rmutil library */ | ||
|
|
||
| static void interp(double x[], double fx[], double *f, double *df) | ||
| { | ||
| int i, j, ni=0; | ||
| double diff1, diff2, tmp1, tmp2, lim1, lim2, tab1[6], tab2[5]; | ||
|
|
||
| tmp1=fabs(x[0]); | ||
| for(i=0;i<5;i++){ | ||
| tmp2=fabs(x[i]); | ||
| if(tmp2<tmp1){ | ||
| ni=i; | ||
| tmp1=tmp2;} | ||
| tab1[i]=fx[i]; | ||
| tab2[i]=fx[i];} | ||
| *f=fx[ni--]; | ||
| for(j=0;j<4;j++){ | ||
| for(i=0;i<=4-j;i++){ | ||
| lim1=x[i]; | ||
| lim2=x[i+j+1]; | ||
| diff1=tab1[i+1]-tab2[i]; | ||
| diff2=lim1-lim2; | ||
| if(diff2==0.0)return; | ||
| diff2=diff1/diff2; | ||
| tab2[i]=lim2*diff2; | ||
| tab1[i]=lim1*diff2;} | ||
| *df=2*ni<(2-j)?tab1[ni+1]:tab2[ni--]; | ||
| *f+=*df;}} | ||
|
|
||
| #define FCN(x) ((*fcn)(x)) | ||
|
|
||
| static double evalfn(double (*fcn)(double), int n) | ||
| { | ||
| double x, nn, tmpsum, pnt1, pnt2; | ||
| static double sum; | ||
| int i,j; | ||
|
|
||
| if (n==1){ | ||
| sum=FCN(0.5); | ||
| return(sum);} | ||
| else { | ||
| for(i=1,j=1;j<n-1;j++) i*=3; | ||
| nn=i; | ||
| pnt1=1/(3.0*nn); | ||
| pnt2=2.0*pnt1; | ||
| x=0.5*pnt1; | ||
| tmpsum=0.0; | ||
| for(j=1;j<=i;j++){ | ||
| tmpsum+=FCN(x); | ||
| x+=pnt2; | ||
| tmpsum+=FCN(x); | ||
| x+=pnt1;} | ||
| sum=(sum+tmpsum/nn)/3.0; | ||
| return(sum);}} | ||
|
|
||
| static double romberg(double (*fcn)(double), double eps) | ||
| { | ||
| int j,j1; | ||
| double sum,errsum=0,x[16],fx[16]; | ||
|
|
||
| x[0]=1.0; | ||
| x[1]=x[0]/9.0; | ||
| x[2]=x[1]/9.0; | ||
| x[3]=x[2]/9.0; | ||
| x[4]=x[3]/9.0; | ||
| x[5]=x[4]/9.0; | ||
| x[6]=x[5]/9.0; | ||
| x[7]=x[6]/9.0; | ||
| x[8]=x[7]/9.0; | ||
| x[9]=x[8]/9.0; | ||
| x[10]=x[9]/9.0; | ||
| x[11]=x[10]/9.0; | ||
| x[12]=x[11]/9.0; | ||
| x[13]=x[12]/9.0; | ||
| x[14]=x[13]/9.0; | ||
| x[15]=x[14]/9.0; | ||
|
|
||
| for(j=0;j<16;j++){ | ||
| j1=j+1; | ||
| fx[j]=evalfn(fcn,j1); | ||
| if(j1>=5){ | ||
| interp(&x[j1-5],&fx[j1-5],&sum,&errsum); | ||
| if(fabs(errsum)<eps*fabs(sum))return(sum);} | ||
| x[j1]=x[j]/9.0; | ||
| fx[j1]=fx[j];} | ||
| sum=R_NaN; | ||
| return(sum);} | ||
|
|
||
| /* density of a stable distribution */ | ||
| void stable(int *n, double *y, double *beta, double *alpha, int *npt, | ||
| double *up, double *eps, int *type, int *err, double *ffy) | ||
| { | ||
| int i, j; | ||
| double h, s, *eta, *seta, *ceta, *sa; | ||
| *err=0; | ||
| eta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| seta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| ceta=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| sa=(double*)R_alloc((size_t)(*n),sizeof(double)); | ||
| nn=*n; | ||
| if(!eta||!seta||!ceta||!sa){ | ||
| *err=1; | ||
| return;} | ||
| for(i=0;i<*n;i++){ | ||
| ffy[i]=0.0; | ||
| eta[i]=beta[i]*(1.0-fabs(1.0-alpha[i]))*M_PI/2.0; | ||
| seta[i]=sin(eta[i]); | ||
| ceta[i]=cos(eta[i]);} | ||
| if(*type==1){ | ||
| *npt=*npt-*npt%2; | ||
| h=*up/ *npt; | ||
| for(j=0;j<=*npt;j++){ | ||
| s=(*npt-j)*h; | ||
| for(i=0;i<*n;i++){ | ||
| sa[i]=pow(s,alpha[i]); | ||
| ffy[i]=ffy[i]+(4-2*(j%2==0)-(j==1||j==*npt))*cos(-y[i]*s+sa[i]*seta[i])*exp(-sa[i]*ceta[i]);}} | ||
| for(i=0;i<*n;i++)ffy[i]=ffy[i]*h/3.0/M_PI;} | ||
| else { | ||
| for(i=0;i<*n;i++){ | ||
| alphai=alpha[i]; | ||
| yi=y[i]; | ||
| setai=seta[i]; | ||
| cetai=ceta[i]; | ||
| ffy[i]=(romberg(fcn1, *eps)+romberg(fcn2, *eps))/M_PI;}}} | ||
|
|
||
| /* cdf of a stable distribution */ | ||
| void pstable_c(int *n, double *y, double *beta, double *alpha, | ||
| double *eps, int *err, double *ffy) | ||
| { | ||
| int i; | ||
| *err=0; | ||
| nn=*n; | ||
| for(i=0;i<*n;i++){ | ||
| ffy[i]=0.0; | ||
| etai=beta[i]*(1.0-fabs(1.0-alpha[i]))*M_PI/2.0; | ||
| setai=sin(etai); | ||
| cetai=cos(etai); | ||
| alphai=alpha[i]; | ||
| yyi=y[i]; | ||
| if(etai==0.&&yyi==0) | ||
| ffy[i]=0.5; | ||
| else ffy[i]=0.5+(romberg(fcn3, *eps)+romberg(fcn4, *eps))/M_PI;}} | ||
|
|
| @@ -0,0 +1,22 @@ | ||
| #include <stdlib.h> // for NULL | ||
| #include <R_ext/Rdynload.h> | ||
|
|
||
| /* FIXME: | ||
| Check these declarations against the C/Fortran source code. | ||
| */ | ||
|
|
||
| /* .C calls */ | ||
| extern void pstable_c(void *, void *, void *, void *, void *, void *, void *); | ||
| extern void stable(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); | ||
|
|
||
| static const R_CMethodDef CEntries[] = { | ||
| {"pstable_c", (DL_FUNC) &pstable_c, 7}, | ||
| {"stable", (DL_FUNC) &stable, 10}, | ||
| {NULL, NULL, 0} | ||
| }; | ||
|
|
||
| void R_init_stable(DllInfo *dll) | ||
| { | ||
| R_registerRoutines(dll, CEntries, NULL, NULL, NULL); | ||
| R_useDynamicSymbols(dll, FALSE); | ||
| } |
