problem in generating positive stable random numbers

"X. Cong" <xcong at stat.rice.edu> writes:

Dear all,

I am trying to use the
rstable(n, alpha, beta, gamma = 1, delta = 0, pm = c(0, 1, 2)))

 function to generate positive stable random numbers. For positive stable
distribution, beta==1 and alpha is in (0,1), which defines random variables
with support (0, infinity). So, I used rstable(100, 0.5, 1) for an example.
I found that this gives me some negative numbers. For example,

rstable(10, 0.5, 1)

 [1]   6.3016252 399.3659030  11.2735789   1.9550625  -0.6762333   1.6810761
 [7]   0.9091360   1.9100991  -0.7593737  24.2788471

Does anybody know why this should happen?

Doesn't sound like it should...

From where did you get rstable()? It isn't part of R itself, and

poking around shows at least three different sources for it
(CircStats, fBasics, stable (Lambert/Lindsey, not on CRAN)). I seem to
recall that the third one also has a distribution function which is
almost, but not quite monotone.

The usual advice of taking package problems to the package maintainer
applies.

O__  ---- Peter Dalgaard             ??ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907

From where did you get rstable()? It isn't part of R itself, and

poking around shows at least three different sources for it
(CircStats, fBasics, stable (Lambert/Lindsey, not on CRAN)). I seem to
recall that the third one also has a distribution function which is
almost, but not quite monotone.

The usual advice of taking package problems to the package maintainer
applies. 

I got rstable() from the fBasics package. Thanks.

Dear all,

I am trying to use the
rstable(n, alpha, beta, gamma = 1, delta = 0, pm = c(0, 1, 2)))

function to generate positive stable random numbers. For positive stable
distribution, beta==1 and alpha is in (0,1), which defines random variables
with support (0, infinity). 

Just a remark - rstable() is from R-package  fBasics ...
1) I think the support for beta=1 and alpha = 1/2 is (-1, infinity), 
isn't it? - Then everything is fine.
2) beta ==1 is a difficult value for numerical computations, try also 
beta = 1-1e-8!
3) Use the program stable.exe from 
http://academic2.american.edu/~jpnolan/stable/stable.html for comparison!
4) Take care in which parametrization you work ...

Best regards
Diethelm


Remark for 3)

  STABLE  3.14.02 (2005/02/28)   Serial number        131
  Copyright 1997-2003 John P. Nolan (jpnolan at american.edu)

  Output file: stable.out                     

  Current tolerance settings:
   -1 debug (F)
    1 relative error for pdf ( 0.1200000000E-13)
    2 relative error for cdf ( 0.1200000000E-13)
    3 relative error for quantiles ( 0.1200000000E-13)
    4 alpha and beta rounding ( 0.1000000000E-01)
    5 x tolerance near zeta ( 0.5000000000E-02)
    6 exponential cutoff (  200.0000000    )
    7 peak/strim location tolerance ( 0.1000000000E-13)
    8 strim tolerance ( 0.1000000000E-50)
    9 minimum alpha ( 0.1000000000    )
   10 minimum xtol ( 0.1000000000E-12)
   11 threshold for quantile search ( 0.1000000000E-09)

   8/30/2005        23:59:34.72  

          Simulation of stable random variables
          n=      10   iseed=             -1   iparam=    0
               alpha      beta     gamma     delta
             0.50000   1.00000   1.00000   0.00000

   39.6516320412757    
  0.247648247004956    
   12.3195079629881    
 -0.146866807337654    
  0.640815989111718    
  0.219559949961403    
   3.05027005351445    
 -0.667188061257615    
   4.02760809108209    
 -0.588547074136516

So, I used rstable(100, 0.5, 1) for an example.
I found that this gives me some negative numbers. For example,

 

rstable(10, 0.5, 1)
   

[1]   6.3016252 399.3659030  11.2735789   1.9550625  -0.6762333   1.6810761
[7]   0.9091360   1.9100991  -0.7593737  24.2788471

Does anybody know why this should happen?

Thanks a lot,
Julie

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html