1,937 NEWS.md
View file

Large diffs are not rendered by default.

784 R/bnlr.r
View file

Large diffs are not rendered by default.

718 R/fitdist.r
View file

Large diffs are not rendered by default.

2,035 R/fmr.r
View file

Large diffs are not rendered by default.

2,501 R/gnlr.r
View file

Large diffs are not rendered by default.

2,406 R/gnlr3.r
View file

Large diffs are not rendered by default.

601 R/nlr.r
View file

Large diffs are not rendered by default.

847 R/nordr.r
View file

Large diffs are not rendered by default.

440 R/ordglm.r
View file

Large diffs are not rendered by default.

372 R/rs.r
View file
@@ -1,152 +1,220 @@
#
# gnlm : A Library of Special Functions for Nonlinear Regression
# Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public Licence as published by
# the Free Software Foundation; either version 2 of the Licence, 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 Licence for more details.
#
# You should have received a copy of the GNU General Public Licence
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
# SYNOPSIS
#
# rs2(y, x1, x2, power=c(1,1), weight=rep(1,length(x1)),
# family=normal, iterlim=20)
# rs3(y, x1, x2, x3, power=c(1,1,1), weight=rep(1,length(x1)),
# family=normal, iterlim=20)
#
# DESCRIPTION
#
# Functions to fit two- and three-covariate power-transformed
# response surface models for glms (Box-Tidwell transformation)

### two covariate model
###
rs2 <- function(y, x1, x2, power=c(1,1), weight=rep(1,length(x1)),
family=normal, iterlim=20){
#
# check estimates and data
#
if(any(c(x1,x2)<0))stop("All covariates must be non-negative")
if(length(power)!=2)
stop("Two estimates of power parameters must be supplied\n")
a <- power[1]
b <- power[2]
#
# iterate, calling glm, to obtain power parameter estimates
#
test <- TRUE
i <- 0
while(test){
xx1 <- x1^a
xx2 <- x2^b
u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,
weight=weight)
z1 <- (u$coef[2]*xx1+2*u$coef[4]*xx1^2+u$coef[6]*xx1*xx2)*
log(ifelse(x1==0,1,x1))
z2 <- (u$coef[3]*xx2+2*u$coef[5]*xx2^2+u$coef[6]*xx1*xx2)*
log(ifelse(x2==0,1,x2))
if(any(is.na(c(z1,z2))))
stop(paste("NAs in calculating estimates:",a,b))
u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2+z1+z2,
family=family,weight=weight)
a <- a+u$coef[6]
b <- b+u$coef[7]
if(any(is.na(c(a,b))))stop(paste("NAs in calculating estimates:",a,b))
i <- i+1
test <- ((u$coef[6]^2>0.00001)||(u$coef[7]^2>0.00001))&&(i<iterlim)}
#
# set up final results
#
z <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,weight=weight)
z$df.residual <- z$df.residual-2
z$aic <- z$aic+4
z$powers <- c(a,b)
z$iterations <- i
class(z) <- c("rs",class(z))
return(z)}

### three covariate model
###
rs3 <- function(y, x1, x2, x3, power=c(1,1,1), weight=rep(1,length(x1)),
family=normal, iterlim=20){
#
# check estimates and data
#
if(any(c(x1,x2,x3)<0))stop("All covariates must be non-negative")
if(length(power)!=3)
stop("Three estimates of power parameters must be supplied\n")
a <- power[1]
b <- power[2]
d <- power[3]
#
# iterate, calling glm, to obtain power parameter estimates
#
test <- TRUE
i <- 0
while(test){
xx1 <- x1^a
xx2 <- x2^b
xx3 <- x3^d
xx12 <- xx1*xx2
xx13 <- xx1*xx3
xx23 <- xx2*xx3
u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
xx12+xx13+xx23,family=family,weight=weight)
z1 <- (u$coef[2]*xx1+2*u$coef[5]*xx1^2+u$coef[8]*xx12+
+u$coef[9]*xx13)*log(ifelse(x1==0,1,x1))
z2 <- (u$coef[3]*xx2+2*u$coef[6]*xx2^2+u$coef[8]*xx12+
u$coef[10]*xx23)*log(ifelse(x2==0,1,x2))
z3 <- (u$coef[4]*xx3+2*u$coef[7]*xx3^2+u$coef[9]*xx13+
u$coef[10]*xx23)*log(ifelse(x3==0,1,x3))
if(any(is.na(c(z1,z2,z3))))
stop(paste("NAs in calculating estimates:",a,b,d))
u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
xx12+xx13+xx23+z1+z2+z3,family=family,weight=weight)
a <- a+u$coef[11]
b <- b+u$coef[12]
d <- d+u$coef[13]
if(any(is.na(c(a,b,d))))stop(paste("NAs in calculating estimates:",a,b,d))
i <- i+1
test <- ((u$coef[11]^2>0.00001)||(u$coef[12]^2>0.00001)||
(u$coef[13]^2>0.00001))&&(i<iterlim)}
#
# set up final results
#
z <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+xx12+xx13+xx23,
family=family,weight=weight)
z$df.residual <- z$df.residual-3
z$aic <- z$aic+6
z$powers <- c(a,b,d)
z$iterations <- i
class(z) <- c("rs",class(z))
return(z)}

### print and summary methods
###
print.rs <- function(z,...){
cat("\nPowered transformed response surface\n\n")
cat("Powers:",z$powers,"\n")
cat("Iterations:",z$iterations,"\n")
print.glm(z,...)}

print.summary.rs <- function(z,...){
cat("\nPowered transformed response surface\n\n")
cat("Powers:",z$powers,"\n")
cat("Iterations:",z$iterations,"\n")
print.summary.glm(z,...)}

summary.rs <- function(z,...){
zz <- summary.glm(z,...)
class(zz) <- c("summary.rs",class(zz))
if(!is.null(z$powers))zz$powers <- z$powers
if(!is.null(z$iterations))zz$iterations <- z$iterations
zz}
#
# gnlm : A Library of Special Functions for Nonlinear Regression
# Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public Licence as published by
# the Free Software Foundation; either version 2 of the Licence, 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 Licence for more details.
#
# You should have received a copy of the GNU General Public Licence
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
# SYNOPSIS
#
# rs2(y, x1, x2, power=c(1,1), weight=rep(1,length(x1)),
# family=gaussian, iterlim=20)
# rs3(y, x1, x2, x3, power=c(1,1,1), weight=rep(1,length(x1)),
# family=gaussian, iterlim=20)
#
# DESCRIPTION
#
# Functions to fit two- and three-covariate power-transformed
# response surface models for glms (Box-Tidwell transformation)

### two covariate model
###


#' Two-factor Box-Tidwell Nonlinear Response Surface Models
#'
#' \code{rs2} fits a two-covariate power-transformed response surface by
#' iterating the function, \code{\link{glm}}.
#'
#'
#' @param y Response variable
#' @param x1 First covariate
#' @param x2 Second covariate
#' @param power Initial estimates of the two power transformations
#' @param weight Weight vector
#' @param family glm family
#' @param iterlim Iteration limit
#' @return A list of class, \code{rs}, is returned containing the model and the
#' power estimates.
#' @author J.K. Lindsey
#' @seealso \code{\link{lm}}, \code{\link{glm}}, \code{\link[gnlm]{gnlr}},
#' \code{\link[gnlm]{gnlr3}}, \code{\link[gnlm]{rs3}}
#' @keywords models
#' @examples
#'
#' x1 <- rep(1:4,5)
#' x2 <- rep(1:5,rep(4,5))
#' y <- rpois(20,1+2*sqrt(x1)+3*log(x2)+4*x1+log(x2)^2+2*sqrt(x1)*log(x2))
#' rs2(y, x1, x2, family=poisson)
#'
#' @export rs2
rs2 <- function(y, x1, x2, power=c(1,1), weight=rep(1,length(x1)),
family=gaussian, iterlim=20){
#
# check estimates and data
#
if(any(c(x1,x2)<0))stop("All covariates must be non-negative")
if(length(power)!=2)
stop("Two estimates of power parameters must be supplied\n")
a <- power[1]
b <- power[2]
#
# iterate, calling glm, to obtain power parameter estimates
#
test <- TRUE
i <- 0
while(test){
xx1 <- x1^a
xx2 <- x2^b
u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,
weights=weight)
z1 <- (u$coef[2]*xx1+2*u$coef[4]*xx1^2+u$coef[6]*xx1*xx2)*
log(ifelse(x1==0,1,x1))
z2 <- (u$coef[3]*xx2+2*u$coef[5]*xx2^2+u$coef[6]*xx1*xx2)*
log(ifelse(x2==0,1,x2))
if(any(is.na(c(z1,z2))))
stop(paste("NAs in calculating estimates:",a,b))
u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2+z1+z2,
family=family,weights=weight)
a <- a+u$coef[6]
b <- b+u$coef[7]
if(any(is.na(c(a,b))))stop(paste("NAs in calculating estimates:",a,b))
i <- i+1
test <- ((u$coef[6]^2>0.00001)||(u$coef[7]^2>0.00001))&&(i<iterlim)}
#
# set up final results
#
z <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,weights=weight)
z$df.residual <- z$df.residual-2
z$aic <- z$aic+4
z$powers <- c(a,b)
z$iterations <- i
class(z) <- c("rs",class(z))
return(z)}

### three covariate model
###


#' Three-factor Box-Tidwell Nonlinear Response Surface Models
#'
#' \code{rs3} fits a three-covariate power-transformed response surface by
#' iterating the function, \code{\link{glm}}.
#'
#'
#' @param y Response variable
#' @param x1 First covariate
#' @param x2 Second covariate
#' @param x3 Third covariate
#' @param power Initial estimates of the three power transformations
#' @param weight Weight vector
#' @param family glm family
#' @param iterlim Iteration limit
#' @return A list of class, \code{rs}, is returned containing the model and the
#' power estimates.
#' @author J.K. Lindsey
#' @seealso \code{\link{lm}}, \code{\link{glm}}, \code{\link[gnlm]{gnlr}},
#' \code{\link[gnlm]{gnlr3}}, \code{\link[gnlm]{rs2}}
#' @keywords models
#' @examples
#'
#' x1 <- rep(1:4,5)
#' x2 <- rep(1:5,rep(4,5))
#' x3 <- c(rep(1:3,6),1,2)
#' #y <- rpois(20,1+2*sqrt(x1)+3*log(x2)+1/x3+4*x1+log(x2)^2+1/x3^2+
#' # 2*sqrt(x1)*log(x2)+sqrt(x1)/x3+log(x2)/x3)
#' y <- c(9,11,14,33,11,19,20,27,22,32,24,24,20,28,25,41,26,31,37,34)
#' rs3(y, x1, x2, x3, family=poisson)
#'
#' @export rs3
rs3 <- function(y, x1, x2, x3, power=c(1,1,1), weight=rep(1,length(x1)),
family=gaussian, iterlim=20){
#
# check estimates and data
#
if(any(c(x1,x2,x3)<0))stop("All covariates must be non-negative")
if(length(power)!=3)
stop("Three estimates of power parameters must be supplied\n")
a <- power[1]
b <- power[2]
d <- power[3]
#
# iterate, calling glm, to obtain power parameter estimates
#
test <- TRUE
i <- 0
while(test){
xx1 <- x1^a
xx2 <- x2^b
xx3 <- x3^d
xx12 <- xx1*xx2
xx13 <- xx1*xx3
xx23 <- xx2*xx3
u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
xx12+xx13+xx23,family=family,weights=weight)
z1 <- (u$coef[2]*xx1+2*u$coef[5]*xx1^2+u$coef[8]*xx12+
+u$coef[9]*xx13)*log(ifelse(x1==0,1,x1))
z2 <- (u$coef[3]*xx2+2*u$coef[6]*xx2^2+u$coef[8]*xx12+
u$coef[10]*xx23)*log(ifelse(x2==0,1,x2))
z3 <- (u$coef[4]*xx3+2*u$coef[7]*xx3^2+u$coef[9]*xx13+
u$coef[10]*xx23)*log(ifelse(x3==0,1,x3))
if(any(is.na(c(z1,z2,z3))))
stop(paste("NAs in calculating estimates:",a,b,d))
u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
xx12+xx13+xx23+z1+z2+z3,family=family,weights=weight)
a <- a+u$coef[11]
b <- b+u$coef[12]
d <- d+u$coef[13]
if(any(is.na(c(a,b,d))))stop(paste("NAs in calculating estimates:",a,b,d))
i <- i+1
test <- ((u$coef[11]^2>0.00001)||(u$coef[12]^2>0.00001)||
(u$coef[13]^2>0.00001))&&(i<iterlim)}
#
# set up final results
#
z <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+xx12+xx13+xx23,
family=family,weights=weight)
z$df.residual <- z$df.residual-3
z$aic <- z$aic+6
z$powers <- c(a,b,d)
z$iterations <- i
class(z) <- c("rs",class(z))
return(z)}

### print and summary methods
###
print.rs <- function(z,...){
cat("\nPowered transformed response surface\n\n")
cat("Powers:",z$powers,"\n")
cat("Iterations:",z$iterations,"\n")
#print.glm(z,...)
class(z) <- "glm"
print(z)}

print.summary.rs <- function(z,...){
cat("\nPowered transformed response surface\n\n")
cat("Powers:",z$powers,"\n")
cat("Iterations:",z$iterations,"\n")
#print.summary.glm(z,...)
class(z) <- "summary.glm"
print(z)}

summary.rs <- function(z,...){
zz <- summary.glm(z,...)
class(zz) <- c("summary.rs",class(zz))
if(!is.null(z$powers))zz$powers <- z$powers
if(!is.null(z$iterations))zz$iterations <- z$iterations
zz}


10 README.md
View file
@@ -0,0 +1,10 @@

<!-- README.md is generated from README.Rmd. Please edit README.Rmd -->
[![Travis-CI Build Status](https://travis-ci.org/swihart/gnlm.svg?branch=master)](https://travis-ci.org/swihart/gnlm) [![CRAN\_Status\_Badge](http://www.r-pkg.org/badges/version/gnlm)](https://cran.r-project.org/package=gnlm) ![downloads](http://cranlogs.r-pkg.org/badges/grand-total/gnlm)

`gnlm` R package
================

This package is intended to be the developmental version to the CRAN version of [Jim Lindsey's gnlm](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/gnlm/compare/jim-lindsey-homepage-version-1.0...master?diff=split&name=master).
68 cran-comments.md
View file
@@ -0,0 +1,68 @@
# gnlm R package
Bruce Swihart
Feb 2019

## Re-Submission 1

* fixed an `Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’` NOTE.

## 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.



# gnlm R package
Bruce Swihart
December 2016

## Test environments
* local OS X install: R version 3.3.2 (2016-10-31)
* ubuntu 12.04 (on travis-ci): R version 3.3.1 (2016-06-21)
* win-builder: R Under development (unstable) (2016-12-05 r71733)

## R CMD check results
There were no ERRORs or WARNINGs.

There was 1 NOTE:

* 1st submission (Package was archived on CRAN, see Miscellaneous)

## Downstream dependencies
There are currently no downstream dependencies for this package.

## Miscellaneous
As per the CRAN Repository Policy Version Revision 3577,

* Explain any change in the maintainer’s email address and if possible send confirmation from the previous address (by a separate email to CRAN@R-project.org) or explain why it is not possible.

This Package was archived on CRAN. I am resurrecting it.
I have Jim Lindsey's permission to be maintainer of his packages on CRAN. Currently he has them in .zip files on his homepage: http://www.commanster.eu/rcode.html. He does not have access to the original email address he used when this package was on CRAN, and thus cannot send a separate confirmation email. However, you can contact him at James Lindsey <jlindsey@gen.unimaas.nl>.

He sent this confirmation message on 2016-11-24:

---------- Forwarded message ----------
From: James Lindsey <jlindsey@gen.unimaas.nl>
Date: Thu, Nov 24, 2016 at 9:26 AM
Subject: transferring maintainer role to Bruce Swihart for 'rmutil', 'repeated', 'gnlm', 'growth', 'event', 'stable'
To: CRAN@r-project.org, ligges@statistik.tu-dortmund.de, Kurt.Hornik@wu.ac.at
Cc: bruce.swihart@gmail.com


Dear CRAN,

Bruce Swihart is taking on the role of maintainer for my R-packages which have been available on my homepage, http://www.commanster.eu/rcode.html. Some of these R-packages were on CRAN but have since been archived.
I know it is CRAN repository policy to email confirmation from the previous maintainer's email address (jlindsey@luc.ac.be), but alas I can't due to LUC changing its name.

Please accept this email as confirmation of the maintainer role changing for R packages 'rmutil', 'repeated', 'gnlm', 'growth', 'event', 'stable'.

Regards,
Jim Lindsey
221 man/bnlr.Rd 100755 → 100644
View file
136 man/fit.dist.Rd 100755 → 100644
View file
369 man/fmr.Rd 100755 → 100644
View file

Large diffs are not rendered by default.

25 man/gnlm.Rd
View file

This file was deleted.

419 man/gnlr.Rd 100755 → 100644
View file

Large diffs are not rendered by default.

383 man/gnlr3.Rd 100755 → 100644
View file

Large diffs are not rendered by default.

193 man/nlr.Rd 100755 → 100644
View file
252 man/nordr.Rd 100755 → 100644
View file
98 man/ordglm.Rd 100755 → 100644
View file
83 man/rs2.Rd 100755 → 100644
View file
92 man/rs3.Rd 100755 → 100644
View file
11 src/Makefile
View file

This file was deleted.

8 src/gnlm_init.c
View file
@@ -0,0 +1,8 @@
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>

void R_init_gnlm(DllInfo *dll)
{
R_registerRoutines(dll, NULL, NULL, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}
372 src/stable.c
View file
@@ -1,186 +1,186 @@
/*
* 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(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=0,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;}}