~ubuntu-branches/ubuntu/vivid/r-cran-stabledist/vivid : revision 6

1

require("stabledist")

2

pPareto <- stabledist:::pPareto

3

4

source(system.file("test-tools.R", package = "Matrix"))

5

#-> identical3(), showProc.time(),...

4

source(system.file("test-tools-1.R", package = "Matrix"), keep.source=interactive())

5

#-> identical3(), showProc.time(),...

6

(doExtras <- stabledist:::doExtras())

6

7

8

options(pstable.debug = FALSE)

8

9

options(pstable.debug = TRUE)# want to see when uniroot() is called

15

16

## a "outer vectorized" version:

16

17

pstabALL <- function(x, alpha, beta, ...)

17

18

sapply(alpha, function(alph)

18

sapply(beta, function(bet)

19

sapply(beta, function(bet)

19

20

pstable(x, alph, bet, ...)))

20

21

22

alph.s <- (1:32)/16 # in (0, 2]

22

23

beta.s <- (-16:16)/16 # in [-1, 1]

23

24

stopifnot(pstabALL( Inf, alph.s, beta.s) == 1,

24

25

pstabALL(-Inf, alph.s, beta.s, log.p=TRUE) == -Inf,

25

pstabALL( 0, alph.s, beta = 0) == 0.5,

26

pstabALL( 0, alph.s, beta = 0) == 0.5,

26

27

TRUE)

27

28

29

pdf("pstab-ex.pdf")

30

28

31

##---- log-scale -------------

29

32

r <- curve(pstable(x, alpha=1.8, beta=.9,

30

lower.tail=FALSE, log.p=TRUE),

31

5, 150, n=500,

32

log="x",type="b", cex=.5)

33

curve(stabledist:::pPareto(x, alpha=1.8, beta=.9,

34

lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)

33

lower.tail=FALSE, log.p=TRUE),

34

5, 150, n=500,

35

log="x",type="b", cex=.5)

36

curve(pPareto(x, alpha=1.8, beta=.9,

37

lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)

35

38

##--> clearly potential for improvement!

36

39

37

40

## the less extreme part - of that:

38

41

r <- curve(pstable(x, alpha=1.8, beta=.9,

39

lower.tail=FALSE, log.p=TRUE),

40

1, 50, n=500, log="x")

41

curve(stabledist:::pPareto(x, alpha=1.8, beta=.9, lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)

42

lower.tail=FALSE, log.p=TRUE),

43

1, 50, n=500, log="x")

44

curve(pPareto(x, alpha=1.8, beta=.9, lower.tail=FALSE, log.p=TRUE), add=TRUE, col=2)

42

45

43

## Check that pstable() is the integral of dstable() --- using simple Simpson's rule

46

## Check that pstable() is the integral of dstable() --- using simple Simpson's rule

44

47

45

48

## in it's composite form:

46

## \int_a^b f(x) \, dx\approx \frac{h}{3}

49

## \int_a^b f(x) dx\approx \frac{h}{3}

47

50

## \bigg[ f(x_0) + 2 \sum_{j=1}^{n/2-1}f(x_{2j}) +

48

## + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) +

49

## + f(x_n) \bigg],

51

## + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) +

52

## + f(x_n) \bigg],

50

53

intSimps <- function(fx, h) {

51

54

stopifnot((n <- length(fx)) %% 2 == 0,

52

55

n >= 4, length(h) == 1, h > 0)

57

60

}

58

61

59

62

chk.pd.stable <- function(alpha, beta, xmin=NA, xmax=NA,

60

n = 256, do.plot=TRUE,

61

comp.tol = 1e-13, eq.tol = 1e-3)

63

n = 256, do.plot=TRUE,

64

comp.tol = 1e-13, eq.tol = 1e-3)

62

65

{

63

66

stopifnot(n >= 20)

64

67

if(is.na(xmin)) xmin <- qstable(0.01, alpha, beta)

69

72

fx <- dstable(x, alpha=alpha, beta=beta, tol= comp.tol)

70

73

Fx <- pstable(x, alpha=alpha, beta=beta, tol=2*comp.tol)

71

74

i.ev <- (i <- seq_along(x))[i %% 2 == 0 & i >= max(n/10, 16)]

72

## integrate from x[1] up to x[i] (where i is even);

75

## integrate from x[1] up to x[i] (where i is even);

73

76

## the exact value will be F(x[i]) - F(x[1]) == Fx[i] - Fx[1]

74

77

Fx. <- vapply(lapply(i.ev, seq_len),

75

function(ii) intSimps(fx[ii], h), 0)

78

function(ii) intSimps(fx[ii], h), 0)

76

79

a.eq <- all.equal(Fx., Fx[i.ev] - Fx[1], tol = eq.tol)

77

80

if(do.plot) {

78

## Show the fit

79

plot(x, Fx - Fx[1], type = "l")

80

lines(x[i.ev], Fx., col=adjustcolor("red", 0.5), lwd=3)

81

op <- par(ask=TRUE) ; on.exit(par(op))

82

## show the "residual", i.e., the relative error

83

plot(x[i.ev], 1- Fx./(Fx[i.ev] - Fx[1]),

84

type = "l", xlim = range(x))

85

abline(h=0, lty=3, lwd = .6)

81

## Show the fit

82

plot(x, Fx - Fx[1], type = "l")

83

lines(x[i.ev], Fx., col=adjustcolor("red", 0.5), lwd=3)

84

op <- par(ask=TRUE) ; on.exit(par(op))

85

## show the "residual", i.e., the relative error

86

plot(x[i.ev], 1- Fx./(Fx[i.ev] - Fx[1]),

87

type = "l", xlim = range(x))

88

abline(h=0, lty=3, lwd = .6)

86

89

}

87

90

88

91

if(!isTRUE(a.eq)) stop(a.eq)

89

92

invisible(list(x=x, f=fx, F=Fx, i. = i.ev, F.appr. = Fx.))

90

93

}

91

94

92

pdf("pstab-ex.pdf")

93

95

op <- par(mfrow=2:1, mar = .1+c(3,3,1,1), mgp=c(1.5, 0.6,0))

94

96

95

97

c1 <- chk.pd.stable(.75, -.5, -1, 1.5, eq.tol = .006)

119

121

## now check convexity/concavity :

120

122

f2 <- diff(diff(px))

121

123

stopifnot(f2[x[-c(1,n)] < 0] > 0,

122

f2[x[-c(1,n)] > 0] < 0)

124

f2[x[-c(1,n)] > 0] < 0)

123

125

## and compare with dstable() ... which actually shows dstable() problem:

124

126

fx. <- dstable(x., alpha=1, beta=.01)

125

127

lines(x., fx., col = 2, lwd=3, lty="5111")

128

130

curve(dstable(x, 1., 0.99), -6, 50, log="y")# "uneven" (x < 0); 50 warnings

129

131

curve(dstable(x, 1.001, 0.95), -10, 30, log="y")# much better

130

132

}

131

c5 <- chk.pd.stable(1., 0.99, -6, 50)# -> uniroot

133

showProc.time() #

134

135

if(doExtras) {

136

c5 <- chk.pd.stable(1., 0.99, -6, 50)# -> uniroot

132

137

c6 <- chk.pd.stable(1.001, 0.95, -10, 30)# -> uniroot; 2nd plot *clearly* shows problem

133

138

with(c5, all.equal(F.appr., F[i.] - F[1], tol = 0)) # .00058 on 64-Lnx

134

139

with(c6, all.equal(F.appr., F[i.] - F[1], tol = 0)) # 2.611e-5 on 64-Lnx

135

140

136

141

## right tail:

137

142

try(## FIXME:

138

c1.0 <- chk.pd.stable(1., 0.8, -6, 500)# uniroot; rel.difference = .030

143

c1.0 <- chk.pd.stable(1., 0.8, -6, 500)# uniroot; rel.difference = .030

139

144

)

140

145

## show it more clearly

141

146

curve(pstable(x, alpha=1, beta=0.5), 20, 800, log="x", ylim=c(.97, 1))

148

153

## zoom

149

154

curve(pstable(x, alpha=1.001, beta=0.5), 100, 200, log="x")

150

155

curve(pPareto(x, alpha=1.001, beta=0.5), add=TRUE, col=2, lty=2)

151

## but alpha = 1 is only somewhat better as approximation:

156

## but alpha = 1 is only somewhat better as approximation:

152

157

curve(pstable(x, alpha=1 , beta=0.5), add=TRUE, col=3,

153

158

lwd=3, lty="5131")

154

155

159

showProc.time() #

156

157

c7 <- chk.pd.stable(1.2, -0.2, -40, 30)

160

}

161

162

163

c7 <- chk.pd.stable(1.2, -0.2, -40, 30)

158

164

c8 <- chk.pd.stable(1.5, -0.999, -40, 30)# two warnings

159

165

160

166

showProc.time() #

167

168

### Newly found -- Marius Hofert, 18.Sept. (via qstable):

169

stopifnot(all.equal(qstable(0.6, alpha = 0.5, beta = 1,

170

tol=1e-15, integ.tol=1e-15),

171

2.636426573120147))

172

##--> which can be traced to the first of

173

stopifnot(pstable(q= -1.1, alpha=0.5, beta=1) == 0,

174

pstable(q= -2.1, alpha=0.6, beta=1) == 0)

175

176

## Stable(alpha = 1/2, beta = 1, gamma, delta, pm = 1) <===> Levy(delta, gamma)

177

source(system.file("xtraR", "Levy.R", package = "stabledist"), keep.source=interactive())

178

##-> dLevy(x, mu, c, log) and

179

##-> pLevy(x, mu, c, log.p, lower.tail)

180

181

set.seed(101)

182

show.Acc <- (interactive() && require("Rmpfr"))

183

if(show.Acc) { ## want to see accuracies, do not stop "quickly"

184

format.relErr <- function(tt, cc)

185

format(as.numeric(relErr(tt, cc)), digits = 4, width = 8)

186

}

187

## FIXME: Look why pstable() is so much less accurate than dstable()

188

## even though the integration in dstable() is more delicate in general

189

190

pTOL <- 1e-6 # typically see relErr of 5e-7

191

dTOL <- 1e-14 # typically see relErr of 1.3...3.9 e-15

192

showProc.time()

193

194

## Note that dstable() is more costly than pstable()

195

for(ii in 1:(if(doExtras) 32 else 8)) {

196

Z <- rnorm(2)

197

mu <- sign(Z[1])*exp(Z[1])

198

sc <- exp(Z[2])

199

x <- seq(mu, mu+ sc* 100*rchisq(1, df=3),

200

length.out= if(doExtras) 512 else 32)

201

## dLevy() and pLevy() using only pnorm() are "mpfr-aware":

202

x. <- if(show.Acc) mpfr(x, 256) else x

203

pL <- pLevy(x., mu, sc)

204

pS <- pstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu,

205

pm = 1)

206

dL <- dLevy(x., mu, sc)

207

dS <- dstable(x, alpha=1/2, beta=1, gamma=sc, delta=mu,

208

pm = 1)

209

if(show.Acc) {

210

cat("p: ", format.relErr(pL, pS), "\t")

211

cat("d: ", format.relErr(dL, dS), "\n")

212

} else {

213

cat(ii %% 10)

214

}

215

stopifnot(all.equal(pL, pS, tol = pTOL),

216

all.equal(dL, dS, tol = dTOL))

217

}; cat("\n")

218

showProc.time()## ~ 75 sec (doExtras=TRUE) on lynne (2012-09)