~ubuntu-branches/ubuntu/utopic/r-cran-stabledist/utopic : revision 1.1.2

45

##' @param alpha

46

##' @return

47

.om <- function(gamma,alpha)

48

if(alpha == 1) (2/pi)*log(gamma) else tan(pi*alpha/2)

48

if(alpha == 1) (2/pi)*log(gamma) else tan(pi2*alpha)

49

50

##' @title C_alpha - the tail constant

51

##' @param alpha numeric vector of stable tail parameters, in [0,2]

64

r

65

}

66

67

##' @title tan(pi/2*x), for x in [-1,1] with correct limits

68

##' i.e. tanpi2(-/+ 1) == -/+ Inf

69

##' @param x numeric vector

70

##' @return numeric vector of values tan(pi/2*x)

71

##' @author Martin Maechler

72

tanpi2 <- function(x) {

73

r <- x

74

if(any(i <- x & x == round(x)))# excluding 0

75

r[i] <- (2 - (x[i] %% 4))*Inf

76

io <- which(!i)

77

r[io] <- tan(pi2* x[io])

78

r

79

}

80

81

##' @title cos(pi/2*x), for x in [-1,1] with correct limits

82

##' i.e. cospi2(+- 1) == 0

83

##' @param x numeric vector

84

##' @return numeric vector of values cos(pi/2*x)

85

##' @author Martin Maechler

86

cospi2 <- function(x) {

87

r <- x

88

if(any(i <- x == round(x)))

89

r[i] <- as.numeric(x[i] == 0)# 1 or 0 - iff x \in [-1,1] !

90

io <- which(!i)

91

r[io] <- cos(pi2* x[io])

92

r

93

}

94

67

95

##' According to Nolan's "tail.pdf" paper, where he takes *derivatives*

68

96

##' of the tail approximation 1-F(x) ~ (1+b) C_a x^{-a} to prove

69

97

##' that f(x) ~ a(1+b) C_a x^{-(1+a)} ...

158

186

159

187

## General Case

160

188

if (alpha != 1) { ## 0 < alpha < 2 & |beta| <= 1 from above

161

tanpa2 <- tan(pi*alpha/2)

189

tanpa2 <- tan(pi2*alpha)

162

190

zeta <- -beta * tanpa2

163

191

theta0 <- min(max(-pi2, atan(beta * tanpa2) / alpha), pi2)

164

192

213

241

r

214

242

}

215

243

244

.e.plus <- function(x, eps) x + eps* abs(x)

245

.e.minus<- function(x, eps) x - eps* abs(x)

246

pi2.. <- function(eps) pi2 * (1 - eps) ## == .e.minus(pi/2, eps), slight more efficiently

247

248

216

249

.fct1 <- function(x, zeta, alpha, theta0, tol, subdivisions, zeta.tol,

217

250

verbose = getOption("dstable.debug", default=FALSE))

218

251

{

235

268

## zeta <- -beta * tan(pi*alpha/2)

236

269

## theta0 <- (1/alpha) * atan( beta * tan(pi*alpha/2))

237

270

## x.m.zet <- abs(x - zeta)

271

##-------->>> identically as in .FCT1() for pstable() below: <<<-----------

238

272

a_1 <- alpha - 1

239

273

cat0 <- cos(at0 <- alpha*theta0)

240

274

##' g() is strictly monotone -- Nolan(1997) ["3. Numerical Considerations"]

275

##' alpha >= 1 <==> g() is falling, ie. from Inf --> 0; otherwise growing from 0 to +Inf

241

276

g <- function(th) {

242

277

r <- th

243

278

## g(-pi/2) or g(pi/2) could become NaN --> work around

244

i.bnd <- abs(pi/2 -sign(a_1)*th) < 64*.Machine$double.eps

279

i.bnd <- abs(pi2 -sign(a_1)*th) < 64*.Machine$double.eps

245

280

r[i.bnd] <- 0

246

281

th <- th[io <- !i.bnd]

247

282

att <- at0 + alpha*th ## = alpha*(theta0 + theta)

263

298

## or dstable(1.205, 0.75, -0.5)

264

299

## the 2nd integral was "completely wrong" (basically zero, instead of ..e-5)

265

300

266

## NB: g() is monotone, typically from 0 to +Inf

267

## alpha >= 1 <==> g() is falling, ie. from Inf --> 0; otherwise growing from 0 to +Inf

301

## NB: g() is monotone, see above

268

302

if((alpha >= 1 &&

269

303

((!is.na(g. <- g( pi2 )) && g. > .large.exp.arg) || identical(g(-theta0), 0))) ||

270

304

(alpha < 1 &&

273

307

## ===> g() * exp(-g()) is 0 everywhere

274

308

return(0)

275

309

276

277

th0.. <- function(th0, eps) -th0+ eps* abs(th0)

278

pi2.. <- function(eps) pi2 * (1 - eps)

279

280

g. <- if(alpha >= 1) g(th0..(theta0, 1e-6)) else g(pi2..(1e-6))

310

g. <- if(alpha >= 1) g(.e.plus(-theta0, 1e-6)) else g(pi2..(1e-6))

281

311

if(is.na(g.))# g() is not usable --- FIXME rather use *asymptotic dPareto()?

282

312

if(max(x.m.zet, x.m.zet / abs(x)) < .01)

283

313

return(f.zeta())

301

331

g1.th2 <- g1( theta2 <- pi2..(1e-6) )

302

332

else if((alpha < 1 && g(-theta0) > 1) ||

303

333

(alpha >= 1 && g(-theta0) < 1))

304

g1.th2 <- g1( theta2 <- th0..(theta0, 1e-6) )

334

g1.th2 <- g1( theta2 <- .e.plus(-theta0, 1e-6) )

305

335

else {

306

336

## when alpha ~=< 1 (0.998 e.g.), g(x) is == 0 (numerically) on a wide range;

307

337

## uniroot is not good enough, and we should *increase* -theta0

367

397

##' @param beta 0 < |beta| <= 1

368

398

##' @param tol

369

399

##' @param subdivisions

370

.fct2 <- function(x, beta, tol, subdivisions)

400

.fct2 <- function(x, beta, tol, subdivisions,

401

verbose = getOption("dstable.debug", default=FALSE))

371

402

{

372

403

i2b <- 1/(2*beta)

373

404

p2b <- pi*i2b # = pi/(2 beta)

374

if(abs(ea <- -p2b*x) > .large.exp.arg) ## ==> g(.) == 0 or == Inf

375

return(0)

376

g.x <- exp(ea)

405

ea <- -p2b*x

406

if(is.infinite(ea)) return(0)

407

377

408

##' g() is strictly monotone;

378

##' for beta > 0: increasing from g(-pi/2) = 0 to g(+pi/2) = Inf

379

##' for beta < 0: decreasing from g(-pi/2) = Inf to g(+pi/2) = 0

380

g <- function(th) {

381

h <- p2b+ th # == g'/beta where g' := pi/2 + beta*th

382

g.x * (h/p2b) * exp(h*tan(th))/cos(th)

383

}

384

t0 <- -sign(beta)*pi2# g(t0) == 0 mathematically, but not always numerically

385

if(is.na(g0 <- g(t0))) { ## numerical problem ==> need more careful version of g()

386

## find a left threshold instead of -pi/2 where g(.) is *not* NA

387

lt <- t0

388

f <- 1 - 1/256

389

while(is.na(g(lt))) lt <- lt*f

390

gg <- g

391

g <- function(th) {

392

r <- th

393

r[i <- (beta*(lt - th) > 0)] <- 0

394

r[!i] <- gg(th[!i])

395

r

396

}

409

##' g(u) := original_g(u*pi/2)

410

##' for beta > 0: increasing from g(-1) = 0 to g(+1) = Inf

411

##' for beta < 0: decreasing from g(-1) = Inf to g(+1) = 0

412

##t0 <- -sign(beta)*pi2# g(t0) == 0 mathematically, but not always numerically

413

u0 <- -sign(beta)# g(u0) == 0 mathematically, but not always numerically

414

g <- function(u) {

415

r <- u

416

r[i <- abs(u-u0) < 1e-10] <- 0

417

u <- u[!i]

418

th <- u*pi2

419

h <- p2b+ th # == g'/beta where g' := pi/2 + beta*th = pi/2* (1 + beta*u)

420

r[!i] <- (h/p2b) * exp(ea + h*tanpi2(u)) / cospi2(u)

421

r

397

422

}

398

423

399

## Function to Integrate; th is a non-sorted vector!

400

g2 <- function(th) {

424

## Function to Integrate; u is a non-sorted vector!

425

g2 <- function(u) {

401

426

## g2 = g(.) exp(-g(.))

402

x.exp.m.x( g(th) )

427

x.exp.m.x( g(u) )

403

428

}

404

429

405

430

## We know that the maximum of g2(.) is = exp(-1) = 0.3679 "at" g(.) == 1

406

431

## find that by uniroot :

407

ur <- uniroot(function(th) g(th) - 1, lower = -pi2, upper = pi2, tol = tol)

408

theta2 <- ur$root

409

410

r1 <- .integrate2(g2, lower = -pi2, upper = theta2,

411

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol)

412

r2 <- .integrate2(g2, lower = theta2, upper = pi2,

413

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol)

414

abs(i2b)*(r1 + r2)

415

}

432

ur <- uniroot(function(u) g(u) - 1, lower = -1, upper = 1, tol = tol)

433

u2 <- ur$root

434

435

r1 <- .integrate2(g2, lower = -1, upper = u2,

436

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol)

437

r2 <- .integrate2(g2, lower = u2, upper = 1,

438

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol)

439

440

cc <- pi2*abs(i2b)

441

if(verbose)

442

cat(sprintf(".fct2(%.11g, %.6g,..): c*sum(r1+r2)= %.11g*(%6.4g + %6.4g)= %g\n",

443

x,beta, cc, r1, r2, cc*(r1+r2)))

444

445

cc*(r1 + r2)

446

}## {.fct2}

416

447

417

448

### ------------------------------------------------------------------------------

418

449

453

484

pcauchy(x, lower.tail=lower.tail, log.p=log.p)

454

485

} else {

455

486

456

retValue <- function(F, useLower) {

487

retValue <- function(F, useLower) { ## (vectorized in F)

457

488

if(useLower) {

458

489

if(log.p) log(F) else F

459

490

} else { ## upper: 1 - F

462

493

}

463

494

## General Case

464

495

if (alpha != 1) { ## 0 < alpha < 2 & |beta| <= 1 from above

465

tanpa2 <- tan(pi*alpha/2)

496

tanpa2 <- tan(pi2*alpha)

466

497

zeta <- -beta * tanpa2

467

498

theta0 <- min(max(-pi2, atan(beta * tanpa2) / alpha), pi2)

468

499

473

504

r <- if(lower.tail) (1/2- theta0/pi) else 1/2+ theta0/pi

474

505

if(log.p) log(r) else r

475

506

} else {

507

useLower <-

508

((z > zeta && lower.tail) ||

509

(z < zeta && !lower.tail))

476

510

## FIXME: for alpha > 1 -- the following computes F1 = 1 -c3*r(x)

477

.F1 <- .FCT1(z, zeta, alpha=alpha,

478

theta0= sign(z - zeta)* theta0,

511

## and suffers from cancellation when 1-F1 is used below:

512

giveI <- !useLower && alpha > 1 # if TRUE, .FCT1() return 1-F

513

.F1 <- .FCT1(z, zeta, alpha=alpha, theta0=theta0,

514

giveI = giveI,

479

515

tol = tol, subdivisions = subdivisions)

480

retValue(.F1, useLower =

481

((z > zeta && lower.tail) ||

482

(z < zeta && !lower.tail)))

516

if(giveI)

517

if(log.p) log(.F1) else .F1

518

else retValue(.F1, useLower=useLower)

483

519

}

484

520

}, 0.)

485

521

}

486

522

## Special Case alpha == 1 and -1 <= beta <= 1 (but not = 0) :

487

523

else { ## (alpha == 1) and 0 < |beta| <= 1 from above

524

useL <-

525

if(beta >= 0)

526

lower.tail

527

else {

528

beta <- -beta

529

x <- -x

530

!lower.tail

531

}

532

if(giveI <- !useL && !log.p)

533

useL <- TRUE

488

534

## Loop over all x values:

489

vapply(x, function(z) {

490

if (beta >= 0) {

491

retValue(.FCT2( z, beta = beta,

492

tol = tol, subdivisions = subdivisions),

493

lower.tail)

494

} else {

495

retValue(.FCT2(-z, beta = -beta,

496

tol = tol, subdivisions = subdivisions),

497

! lower.tail)

498

}

499

}, 0.)

535

retValue(vapply(x, function(z)

536

.FCT2(z, beta = beta, tol=tol, subdivisions=subdivisions,

537

giveI = giveI),

538

0.),

539

useLower = useL)

500

540

}

501

541

}

502

542

}## {pstable}

504

544

## ------------------------------------------------------------------------------

505

545

506

546

##' Auxiliary for pstable() (for alpha != 1)

507

.FCT1 <- function(x, zeta, alpha, theta0, tol, subdivisions)

547

.FCT1 <- function(x, zeta, alpha, theta0, giveI, tol, subdivisions,

548

verbose = getOption("pstable.debug", default=FALSE))

508

549

{

509

550

if(is.infinite(x))

510

return(1)

551

return(if(giveI) 0 else 1)

511

552

x.m.zet <- abs(x - zeta)

512

## identically as in .fct1() for dstable():

553

##-------->>> identically as in .fct1() for dstable() above: <<<-----------

554

if(x < zeta) theta0 <- -theta0

555

513

556

a_1 <- alpha - 1

514

557

cat0 <- cos(at0 <- alpha*theta0)

515

## Nolan(1997) shows that g() is montone

558

516

559

g <- function(th) {

517

560

r <- th

518

561

## g(-pi/2) or g(pi/2) could become NaN --> work around

519

i.bnd <- abs(pi/2 -sign(a_1)*th) < 64*.Machine$double.eps

562

i.bnd <- abs(pi2 -sign(a_1)*th) < 64*.Machine$double.eps

520

563

r[i.bnd] <- 0

521

564

th <- th[io <- !i.bnd]

522

565

att <- at0 + alpha*th ## = alpha*(theta0 + theta)

524

567

r

525

568

}

526

569

527

## as g() is montone, the integrand exp(-g(.)) is too -- so the maximum is at the boundary

570

if(verbose) cat(sprintf(".FCT1(%9g, %10g, th0=%.10g, %s..): ",

571

x,zeta, theta0, if(giveI)"giveI=TRUE," else ""))

572

573

## as g() is montone, the integrand exp(-g(.)) is too ==> maximum is at the boundary

574

## however, integration can be inaccuracte when g(.) quickly jumps from Inf to 0

575

## _BUT_ empirically I find that good values l.th / u.th below are *INDEPENDENT* of x,

576

l.th <- .e.plus(-theta0, 1e-6)

577

if(alpha > 1 && g(l.th) == Inf) {

578

ur <- uniroot(function(t) 1-2*(g(t)==Inf), lower=l.th, upper=pi2,

579

f.lower= -1, f.upper= 1, tol = 1e-8)

580

l.th <- ur$root

581

if(verbose) cat(sprintf(" g(-th0 +1e-6)=Inf: unirt(%d it) -> l.th=%.10g ",

582

ur$iter, l.th))

583

}

584

u.th <- .e.minus(pi2, 1e-6)

585

if(alpha < 1 && g(u.th) == Inf) {

586

ur <- uniroot(function(t) 1-2*(g(t)==Inf), lower=l.th, upper=u.th,

587

f.upper= -1, tol = 1e-8)

588

u.th <- ur$root

589

if(verbose) cat(sprintf(" g(pi/2 -1e-6)=Inf: unirt(%d it) -> u.th=%.10g ",

590

ur$iter, u.th))

591

}

528

592

r <- .integrate2(function(th) exp(-g(th)),

529

lower = -theta0, upper = pi2, subdivisions = subdivisions,

593

lower = l.th, upper = u.th, subdivisions = subdivisions,

530

594

rel.tol = tol, abs.tol = tol)

531

595

532

c1 <- if(alpha < 1) 1/2 - theta0/pi else 1

533

c3 <- sign(1-alpha)/pi

534

## = 1 - |.|*r(x) <==> cancellation iff we eventually want 1 - F() -- FIXME

535

c1 + c3* r

536

}

596

if(verbose) cat(sprintf("--> Int r= %.11g\n", r))

597

if(giveI) { ## { ==> alpha > 1 ==> c1 = 1; c3 = -1/pi}

598

## return (1 - F) = 1 - (1 -1/pi * r) = r/pi :

599

r/pi

600

} else {

601

c1 <- if(alpha < 1) 1/2 - theta0/pi else 1

602

c3 <- sign(1-alpha)/pi

603

## = 1 - |.|*r(x) <==> cancellation iff we eventually want 1 - F() -- FIXME

604

c1 + c3* r

605

}

606

} ## {.FCT1}

537

607

538

608

## ------------------------------------------------------------------------------

539

609

542

612

##' @param beta >= 0 here

543

613

##' @param tol

544

614

##' @param subdivisions

545

.FCT2 <- function(x, beta, tol, subdivisions)

615

.FCT2 <- function(x, beta, tol, subdivisions, giveI = FALSE,

616

verbose = getOption("pstable.debug", default=FALSE))

546

617

{

547

618

i2b <- 1/(2*beta)

548

619

p2b <- pi*i2b # = pi/(2 beta)

549

if((ea <- -p2b*x) < -.large.exp.arg) ## ==> g(.) == 0 ==> G2(.) == 1

550

return(1)

551

if(ea > .large.exp.arg) ## ==> g(.) == Inf ==> G2(.) == 0

552

return(0)

553

g.x <- exp(ea)

620

ea <- -p2b*x

621

if(is.infinite(ea))

622

return(R.D.Lval(if(ea < 0) ## == -Inf ==> g(.) == 0 ==> G2(.) == 1

623

1 else 0, ## == +Inf ==> g(.) == Inf ==> G2(.) == 0

624

lower.tail= !giveI))

625

554

626

##' g() is strictly monotone;

555

##' for beta > 0: increasing from g(-pi/2) = 0 to g(+pi/2) = Inf

556

##' for beta < 0: decreasing from g(-pi/2) = Inf to g(+pi/2) = 0

557

g <- function(th) {

558

h <- p2b+ th # == g'/beta where g' := pi/2 + beta*th

559

g.x * (h/p2b) * exp(h*tan(th))/cos(th)

627

##' g(u) := original_g(u*pi/2)

628

##' for beta > 0: increasing from g(-1) = 0 to g(+1) = Inf

629

##' for beta < 0: decreasing from g(-1) = Inf to g(+1) = 0

630

## original_g :

631

## g <- function(th) {

632

## h <- p2b+ th # == g'/beta where g' := pi/2 + beta*th

633

## (h/p2b) * exp(ea + h*tan(th)) / cos(th)

634

## }

635

##t0 <- -pi2# g(t0) == 0 mathematically, but not always numerically

636

u0 <- -1 # g(u0) == 0 mathematically, but not always numerically

637

g <- function(u) {

638

r <- u

639

r[i <- abs(u-u0) < 1e-10] <- 0

640

u <- u[!i]

641

th <- u*pi2

642

h <- p2b+ th # == g'/beta where g' := pi/2 + beta*th = pi/2* (1 + beta*u)

643

r[!i] <- (h/p2b) * exp(ea + h*tanpi2(u)) / cospi2(u)

644

r

560

645

}

561

t0 <- -sign(beta)*pi2# g(t0) == 0 mathematically, but not always numerically

562

if(is.na(g0 <- g(t0))) { ## numerical problem ==> need more careful version of g()

563

## find a left threshold instead of -pi/2 where g(.) is *not* NA

564

lt <- t0

565

f <- 1 - 1/256

566

while(is.na(g(lt))) lt <- lt*f

567

gg <- g

568

g <- function(th) {

569

r <- th

570

r[i <- (beta*(lt - th) > 0)] <- 0

571

r[!i] <- gg(th[!i])

572

r

573

}

646

647

if(verbose)

648

cat(sprintf(".FCT2(%.11g, %.6g, %s..): ",

649

x,beta, if(giveI) "giveI=TRUE," else ""))

650

651

## g(-u0) == +Inf {at other end}, mathematically ==> exp(-g(.)) == 0

652

## in the outer tails, the numerical integration can be inaccurate,

653

## because g(.) jumps from 0 to Inf, but is 0 almost always

654

## <==> g1(.) = exp(-g(.)) jumps from 1 to 0 and is 1 almost everywhere

655

## ---> the integration "does not see the 0" and returns too large..

656

u. <- 1

657

if(g(uu <- .e.minus(u., 1e-6)) == Inf) {

658

ur <- uniroot(function(t) 1-2*(g(t)==Inf), lower=-1, upper= uu,

659

f.lower= +1, f.upper= -1, tol = 1e-8)

660

u. <- ur$root

661

if(verbose) cat(sprintf(" g(%g)=Inf: unirt(%d it) -> u.=%.10g",

662

uu, ur$iter, u.))

574

663

}

575

664

576

665

##' G2(.) = exp(-g(.)) is strictly monotone .. no need for 'theta2' !

577

.integrate2(function(th) exp(-g(th)), lower = -pi2, upper = pi2,

578

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol) / pi

579

}

666

G2 <- if(giveI) function(u) expm1(-g(u)) else function(u) exp(-g(u))

667

r <- .integrate2(G2, lower = -1, upper = u.,

668

subdivisions = subdivisions, rel.tol = tol, abs.tol = tol) / 2

669

if(verbose) cat(sprintf("--> Int r= %.11g\n", r))

670

if(giveI) -r else r

671

}## {.FCT2}

580

672

581

673

### ------------------------------------------------------------------------------

582

674