~ubuntu-branches/ubuntu/precise/r-cran-stabledist/precise : revision 1

1

2

## This R package is free software; you can redistribute it and/or

3

## modify it under the terms of the GNU Library General Public

4

## License as published by the Free Software Foundation; either

5

## version 2 of the License, or (at your option) any later version.

6

##

7

## This R package is distributed in the hope that it will be useful,

8

## but WITHOUT ANY WARRANTY; without even the implied warranty of

9

## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

10

## GNU Library General Public License for more details.

11

##

12

## You should have received a copy of the GNU Library General

13

## Public License along with this library; if not, write to the

14

## Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,

15

## MA 02111-1307 USA

16

17

18

################################################################################

19

## FUNCTIONS: DESCRIPTION:

20

## dstable Returns density for stable DF

21

## pstable Returns probabilities for stable DF

22

## qstable Returns quantiles for stable DF

23

## rstable Returns random variates for stable DF

24

## UTILITY FUNCTION DESCRIPTION:

25

## .integrate2 Integrates internal functions for *stable

26

################################################################################

27

28

29

##==============================================================================

30

31

### MM TODO:

32

33

## 0) All d, p, q, q -- have identical parts

34

## a) 'parametrization' (pm) check

35

## b) checking, alpha, beta,

36

## c) subdivisions etc {all but rstable}

37

## --- to do: "Fix" these in dstable(), then copy/paste to others

38

39

##==============================================================================

40

41

##' @title omega() according to Lambert & Lindsey (1999), p.412

42

##' @param gamma

43

##' @param alpha

44

##' @return

45

.om <- function(gamma,alpha)

46

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

47

48

dstable <- function(x, alpha, beta,

49

gamma = 1, delta = 0, pm = 0, log = FALSE,

50

tol = 16*.Machine$double.eps, subdivisions = 1000)

51

{

52

## Original implemented by Diethelm Wuertz;

53

## Changes for efficiency and accuracy by Martin Maechler

54

55

## Description:

56

## Returns density for stable DF

57

58

## Details:

59

## The function uses the approach of J.P. Nolan for general

60

## stable distributions. Nolan derived expressions in form

61

## of integrals based on the charcteristic function for

62

## standardized stable random variables. These integrals

63

## can be numerically evaluated.

64

65

## Arguments:

66

## alpha = index of stability, in the range (0,2]

67

## beta = skewness, in the range [-1, 1]

68

## gamma = scale, in the range (0, infinity)

69

## delta = location, in the range (-infinity, +infinity)

70

## param = type of parmeterization

71

72

## Note: S+ compatibility no longer considered (explicitly)

73

74

## Parameter Check:

75

## NB: (gamma, delta) can be *vector*s (vectorized along x)

76

stopifnot( 0 < alpha, alpha <= 2, length(alpha) == 1,

77

-1 <= beta, beta <= 1, length(beta) == 1,

78

0 <= gamma, length(pm) == 1, pm %in% 0:2,

79

tol > 0, subdivisions > 0)

80

81

## Parameterizations:

82

if (pm == 1) {

83

delta <- delta + beta*gamma * .om(gamma,alpha)

84

} else if (pm == 2) {

85

delta <- delta - alpha^(-1/alpha)*gamma*stableMode(alpha, beta)

86

gamma <- alpha^(-1/alpha) * gamma

87

} ## else pm == 0

88

89

## Shift and Scale:

90

x <- (x - delta) / gamma

91

ans <-

92

## Special Cases:

93

if (alpha == 2) {

94

dnorm(x, mean = 0, sd = sqrt(2))

95

} else if (alpha == 1 && beta == 0) {

96

dcauchy(x)

97

} else {

98

99

## General Case

100

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

101

tanpa2 <- tan(pi*alpha/2)

102

zeta <- -beta * tanpa2

103

theta0 <- atan(beta * tanpa2) / alpha

104

## Loop over all x values:

105

unlist(lapply(x, function(z)

106

.fct1(z, zeta, alpha=alpha, theta0 = theta0,

107

tol=tol, subdivisions=subdivisions)))

108

}

109

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

110

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

111

## Loop over all x values:

112

unlist(lapply(x, function(z) {

113

if (z >= 0) {

114

.fct2( z , beta, tol=tol, subdivisions=subdivisions)

115

} else {

116

.fct2(-z, -beta, tol=tol, subdivisions=subdivisions)

117

}

118

}))

119

}

120

}

121

122

ans <- ans/gamma

123

## Return:

124

if (log) log(ans) else ans

125

}

126

127

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

128

129

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

130

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

131

{

132

## -- Integrand for dstable() --

133

134

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

135

f.zeta <- function()

136

gamma(1+1/alpha)*cos(theta0) / (pi*(1+zeta^2)^(1/(2*alpha)))

137

138

## Modified: originally was if (z == zeta),

139

## then (D.W.) if (x.m.zet < 2 * .Machine$double.eps)

140

## then (M.M.) if (x.m.zet <= 1e-5 * abs(x))

141

if(x.m.zet < 1e-10)

142

return(f.zeta())

143

## the real check should be about the feasibility of g() below, or its integration

144

145

if(x < zeta) theta0 <- -theta0

146

147

## constants ( independent of integrand g1(th) = g*exp(-g) ):

148

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

149

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

150

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

151

a_1 <- alpha - 1

152

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

153

154

g <- function(th) (cat0 * cos(th) * (x.m.zet/sin(at0+alpha*th))^alpha)^(1/a_1) * cos(at0+ a_1*th)

155

156

## Function to Integrate:

157

g1 <- function(th) ## and (xarg, alpha, beta) => (zeta,at0,g0,.....)

158

{

159

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

160

v <- g(th)

161

## making sure that we don't get Inf * 0 => NaN

162

r <- v * exp(-v)

163

if(any(L <- v > 1000)) ## actually, v > -(.Machine$double.min.exp * log(2)) == 708.396...

164

r[L] <- 0

165

r

166

}

167

168

c2 <- ( alpha / (pi*abs(a_1)*x.m.zet) )

169

## Result = c2 * \int_{-t0}^{pi/2} g1(u) du

170

## where however, g1(.) may look to be (almost) zero almost everywhere and just have a small peak

171

## ==> Split the integral into two parts of two intervals (t0, t_max) + (t_max, pi/2)

172

173

## However, this may still be bad, e.g., for dstable(71.61531, alpha=1.001, beta=0.6),

174

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

175

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

176

## FIXME --- Lindsey uses "Romberg" integration -- maybe we must split even more

177

178

179

g. <- if(alpha >= 1) g(-theta0+ 1e-6*abs(theta0)) else g(pi/2 -1e-6*(pi/2))

180

if(is.na(g.) || identical(g., 0))# g() is not usable

181

return(f.zeta())

182

183

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

184

## find that by uniroot :

185

ur <- uniroot(function(th) g(th) - 1, lower = -theta0, upper = pi/2, tol = tol)

186

theta2 <- ur$root

187

188

## now, because the peak may be extreme, we find two more intermediate values

189

## NB: Max = 0.3679 ==> 1e-4 is a very small fraction of that

190

## to the left:

191

th1 <- uniroot(function(th) g1(th) - 1e-4, lower = -theta0, upper = theta2,

192

tol = tol)$root

193

if((do4 <- g1(pi/2) < 1e-4))

194

## to the right:

195

th3 <- uniroot(function(th) g1(th) - 1e-4, lower = theta2, upper = pi/2,

196

tol = tol)$root

197

198

Int <- function(a,b)

199

.integrate2(g1, lower = a, upper = b,

200

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

201

202

r1 <- Int(-theta0, th1)

203

r2 <- Int( th1, theta2)

204

if(do4) {

205

r3 <- Int( theta2, th3)

206

r4 <- Int( th3, pi/2)

207

} else {

208

r3 <- Int( theta2, pi/2)

209

r4 <- 0

210

}

211

if(verbose)

212

cat(sprintf(".fct1(%11g,%10g,..): c2*sum(r[1..4])= %11g*(%9.4g + %9.4g + %9.4g + %9.4g)= %g\n",

213

x,zeta, c2, r1,r2,r3,r4, c2*(r1+r2+r3+r4)))

214

c2*(r1+r2+r3+r4)

215

}

216

217

218

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

219

220

##' is only used when alpha == 1 (!)

221

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

222

{

223

## Integration:

224

225

i2b <- 1/(2*beta)

226

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

227

pi2 <- pi/2

228

g. <- exp(-p2b*xarg)

229

230

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

231

g2 <- function(x) { ## and (xarg, beta)

232

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

233

v <- g / (p2b*cos(x)) * exp(g*tan(x))

234

g <- g. * v

235

gval <- g * exp(-g)

236

if(any(ina <- is.na(gval))) gval[ina] <- 0 ## replace NA at pi/2

237

gval

238

}

239

240

theta2 <- optimize(g2, lower = -pi2, upper = pi2,

241

maximum = TRUE, tol = tol)$maximum

242

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

243

subdivisions = subdivisions,

244

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

245

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

246

subdivisions = subdivisions,

247

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

248

abs(i2b)*(r1 + r2)

249

}

250

251

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

252

253

254

pstable <- function(q, alpha, beta, gamma = 1, delta = 0, pm = 0,

255

tol = 16*.Machine$double.eps, subdivisions = 1000)

256

{

257

## A function implemented by Diethelm Wuertz

258

259

## Description:

260

## Returns probability for stable DF

261

262

x <- q

263

## Parameter Check:

264

## NB: (gamma, delta) can be *vector*s (vectorized along x)

265

stopifnot( 0 < alpha, alpha <= 2, length(alpha) == 1,

266

-1 <= beta, beta <= 1, length(beta) == 1,

267

0 <= gamma, length(pm) == 1, pm %in% 0:2,

268

tol > 0, subdivisions > 0)

269

270

## Parameterizations:

271

if (pm == 1) {

272

delta <- delta + beta*gamma * .om(gamma,alpha)

273

} else if (pm == 2) {

274

delta <- delta - alpha^(-1/alpha)*gamma*stableMode(alpha, beta)

275

gamma <- alpha^(-1/alpha) * gamma

276

} ## else pm == 0

277

278

## Shift and Scale:

279

x <- (x - delta) / gamma

280

281

## Return directly

282

## ------ first, special cases:

283

if (alpha == 2) {

284

pnorm(x, mean = 0, sd = sqrt(2))

285

} else if (alpha == 1 && beta == 0) {

286

pcauchy(x)

287

} else {

288

289

## General Case

290

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

291

tanpa2 <- tan(pi*alpha/2)

292

zeta <- -beta * tanpa2

293

theta0 <- atan(beta * tanpa2) / alpha

294

295

## Loop over all x values:

296

unlist(lapply(x, function(z) {

297

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

298

299

if (z.m.zet < 2 * .Machine$double.eps)

300

## FIXME? same problem as dstable

301

(1/2- theta0/pi)

302

else if (z > zeta)

303

.FCT1(z.m.zet, alpha=alpha, theta0= theta0,

304

tol = tol, subdivisions = subdivisions)

305

else ## (z < zeta)

306

1 - .FCT1(z.m.zet, alpha=alpha, theta0= -theta0,

307

tol = tol, subdivisions = subdivisions)

308

}))

309

}

310

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

311

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

312

## Loop over all x values:

313

unlist(lapply(x, function(z) {

314

if (beta >= 0) {

315

.FCT2(xarg = z, beta = beta,

316

tol = tol, subdivisions = subdivisions)

317

} else {

318

1 - .FCT2(xarg = -z, beta = -beta,

319

tol = tol, subdivisions = subdivisions)

320

}

321

}))

322

}

323

}

324

}

325

326

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

327

328

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

329

{

330

a_1 <- alpha - 1

331

aa1 <- alpha/a_1

332

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

333

g0 <- x.m.zeta^aa1

334

335

## Function to integrate:

336

G1 <- function(x) { ## (xarg, alpha, beta)

337

v <- (cat0*cos(x))^(1/a_1) * sin(at0+ alpha*x)^-aa1 * cos(at0+a_1*x)

338

exp(-(g0 * v))

339

}

340

341

theta2 <- optimize(G1, lower = -theta0, upper = pi/2,

342

maximum = TRUE, tol = tol)$maximum

343

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

344

c3 <- sign(1-alpha)/pi

345

r1 <- .integrate2(G1, lower = -theta0,

346

upper = theta2, subdivisions = subdivisions,

347

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

348

r2 <- .integrate2(G1, lower = theta2,

349

upper = pi/2, subdivisions = subdivisions,

350

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

351

c1 + c3*(r1+r2)

352

}

353

354

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

355

356

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

357

{

358

i2b <- 1/(2*beta)

359

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

360

pi2 <- pi/2

361

g. <- exp(-p2b*xarg)

362

363

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

364

G2 <- function(x) { ## and (xarg, beta)

365

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

366

v <- g / (p2b*cos(x)) * exp(g*tan(x))

367

g <- g. * v

368

gval <- exp(-g)

369

if(any(ina <- is.na(gval))) gval[ina] <- 0 ## replace NA at pi/2

370

gval

371

}

372

373

## Integration:

374

theta2 <- optimize(G2, lower = -pi2, upper = pi2,

375

maximum = TRUE, tol = tol)$maximum

376

r1 <- .integrate2(G2, lower = -pi2,

377

upper = theta2, subdivisions = subdivisions,

378

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

379

r2 <- .integrate2(G2, lower = theta2,

380

upper = pi2, subdivisions = subdivisions,

381

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

382

(r1+r2)/pi

383

}

384

385

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

386

387

388

qstable <- function(p, alpha, beta, gamma = 1, delta = 0, pm = 0,

389

tol = .Machine$double.eps^0.25, maxiter = 1000,

390

integ.tol = 1e-7, subdivisions = 200)

391

{

392

## A function implemented by Diethelm Wuertz

393

394

## Description:

395

## Returns quantiles for stable DF

396

397

## Parameter Check:

398

## NB: (gamma, delta) can be *vector*s (vectorized along x)

399

stopifnot( 0 < alpha, alpha <= 2, length(alpha) == 1,

400

-1 <= beta, beta <= 1, length(beta) == 1,

401

0 <= gamma, length(pm) == 1, pm %in% 0:2,

402

tol > 0, subdivisions > 0)

403

404

## Parameterizations:

405

if (pm == 1) {

406

delta <- delta + beta*gamma * .om(gamma,alpha)

407

} else if (pm == 2) {

408

delta <- delta - alpha^(-1/alpha)*gamma*stableMode(alpha, beta)

409

gamma <- alpha^(-1/alpha) * gamma

410

} ## else pm == 0

411

412

result <-

413

## Special Cases:

414

if (alpha == 2)

415

qnorm(p, mean = 0, sd = sqrt(2))

416

else if (alpha == 1 && beta == 0)

417

qcauchy(p)

418

else if (abs(alpha-1) < 1) { ## -------------- 0 < alpha < 2 ---------------

419

.froot <- function(x, p) {

420

pstable(q = x, alpha=alpha, beta=beta, pm = 0,

421

tol=integ.tol, subdivisions=subdivisions) - p

422

}

423

## Calculate:

424

unlist(lapply(p, function(pp) {

425

if (beta < 0) {

426

xmin <- -(1-pp)/pp

427

## xmax = pp/(1-pp)

428

if (pp < 0.5) {

429

xmax <- qnorm(pp, mean = 0, sd = sqrt(2))

430

} else {

431

xmax <- qcauchy(pp)

432

}

433

}

434

if (beta > 0 ) {

435

## xmin = -(1-pp)/pp

436

if (pp < 0.5) {

437

xmin <- qcauchy(pp)

438

} else {

439

xmin <- qnorm(pp, mean = 0, sd = sqrt(2))

440

}

441

xmax <- pp/(1-pp)

442

}

443

if (beta == 0 ) {

444

## xmin = -(1-pp)/pp

445

if (pp < 0.5) {

446

xmin <- qcauchy(pp)

447

} else {

448

xmin <- qnorm(pp, mean = 0, sd = sqrt(2))

449

}

450

## xmax = pp/(1-pp)

451

if (pp < 0.5) {

452

xmax <- qnorm(pp, mean = 0, sd = sqrt(2))

453

} else {

454

xmax <- qcauchy(pp)

455

}

456

}

457

root <- NA

458

counter <- 0

459

while (is.na(root)) {

460

root <- .unirootNA(.froot, interval = c(xmin, xmax), p = pp,

461

tol=tol, maxiter=maxiter)

462

counter <- counter + 1

463

xmin <- xmin-2^counter

464

xmax <- xmax+2^counter

465

}

466

root

467

}))

468

}

469

470

## Result:

471

result * gamma + delta

472

}

473

474

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

475

476

477

rstable <- function(n, alpha, beta, gamma = 1, delta = 0, pm = 0)

478

{

479

## Description:

480

## Returns random variates for stable DF

481

482

## slightly amended along nacopula::rstable1

483

484

## Parameter Check:

485

## NB: (gamma, delta) can be *vector*s (vectorized along x)

486

stopifnot( 0 < alpha, alpha <= 2, length(alpha) == 1,

487

-1 <= beta, beta <= 1, length(beta) == 1,

488

0 <= gamma, length(pm) == 1, pm %in% 0:2)

489

490

## Parameterizations:

491

if (pm == 1) {

492

delta <- delta + beta*gamma * .om(gamma,alpha)

493

} else if (pm == 2) {

494

delta <- delta - alpha^(-1/alpha)*gamma*stableMode(alpha, beta)

495

gamma <- alpha^(-1/alpha) * gamma

496

} ## else pm == 0

497

498

## Calculate uniform and exponential distributed random numbers:

499

theta <- pi * (runif(n)-1/2)

500

w <- -log(runif(n))

501

502

result <-

503

## If alpha is equal 1 then:

504

if (alpha == 1 & beta == 0) {

505

rcauchy(n)

506

## Otherwise, if alpha is different from 1:

507

} else {

508

## FIXME: learn from nacopula::rstable1R()

509

b.tan.pa <- beta*tan(pi*alpha/2)

510

theta0 <- atan(b.tan.pa) / alpha ## == \theta_0

511

c <- (1+b.tan.pa^2)^(1/(2*alpha))

512

a.tht <- alpha*(theta+theta0)

513

r <- ( c*sin(a.tht)/

514

(cos(theta))^(1/alpha) ) *

515

(cos(theta-a.tht)/w)^((1-alpha)/alpha)

516

## Use Parametrization 0:

517

r - b.tan.pa

518

}

519

520

## Result:

521

result * gamma + delta

522

}

523

524

525

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

526

527

528

##' Numerically Integrate -- basically the same as R's integrate()

529

##' --------------------- main difference: no errors, but warnings

530

.integrate2 <- function(f, lower, upper, subdivisions, rel.tol, abs.tol, ...)

531

{

532

## Originally implemented by Diethelm Wuertz -- without *any* warnings

533

534

if (class(version) != "Sversion") {

535

## R:

536

f <- match.fun(f)

537

ff <- function(x) f(x, ...)

538

wk <- .External("call_dqags", ff,

539

rho = environment(), as.double(lower),

540

as.double(upper), as.double(abs.tol),

541

as.double(rel.tol), limit = as.integer(subdivisions),

542

PACKAGE = "base")[c("value","ierr")]

543

iErr <- wk[["ierr"]]

544

if(iErr == 6) stop("the input is invalid")

545

if(iErr > 0)

546

## NB: "roundoff error ..." happens many times

547

warning(switch(iErr + 1, "OK",

548

"maximum number of subdivisions reached",

549

"roundoff error was detected",

550

"extremely bad integrand behaviour",

551

"roundoff error is detected in the extrapolation table",

552

"the integral is probably divergent"))

553

wk[[1]]

554

555

} else {

556

## SPlus:

557

integrate(f, lower, upper, subdivisions, rel.tol, abs.tol, ...)[[1]]

558

}

559

}

560

561

562

################################################################################

563