1

I'm trying to generate random values from a given distribution given by the following function in R with the inverse cdf method found here described after the function

sigma=function(x){
  
  rc=0.5
  
  y=(x/rc) ; 
  
  C=(1/rc) ;
  
  Fun = log(1+C)-( C/(1+C) ) ; 
  
##### terms  
  
  a = 1/( y^2 - 1) ; b = 2/( sqrt(y^2 + 1 ) ) ;  c = atan( sqrt((y - 1)/y+1) ) ; 
  
  ff = a*(1 - b*c) ; 
  
##### My function   
  
  cst=2*pi
  
  sigma1 = (ff)/(cst*(rc^2)*Fun) ## function to generate random values
  
  # return(sigma1)
  
}

######## inverse CDF method from "found here" link

den<-sigma

#calculates the cdf by numerical integration
cdf<-function(x) integrate(den,-Inf,x)[[1]]

#inverts the cdf
inverse.cdf<-function(x,cdf,starting.value=0){
 lower.found<-FALSE
 lower<-starting.value
 while(!lower.found){
  if(cdf(lower)>=(x-.000001))
   lower<-lower-(lower-starting.value)^2-1
  else
   lower.found<-TRUE
 }
 upper.found<-FALSE
 upper<-starting.value
 while(!upper.found){
  if(cdf(upper)<=(x+.000001))
   upper<-upper+(upper-starting.value)^2+1
  else
   upper.found<-TRUE
 }
 uniroot(function(y) cdf(y)-x,c(lower,upper))$root
}

#generates 1000 random variables of distribution 'den'
vars<-apply(matrix(runif(1000)),1,function(x) inverse.cdf(x,cdf))
hist(vars)

And then the following error appears:

Error in integrate(den, -Inf, x) : the integral is probably divergent
Called from: integrate(den, -Inf, x)
``

Unfortunately I can't understand where the problem is in my function. However I tried another approach as follows: ([which can be found here][2])

    library(distr)
    p <- sigma  # probability density function
    dist <-AbscontDistribution(d=p)  # signature for a dist with pdf ~ p
    rdist <- r(dist)                 # function to create random variates from p
    
    set.seed(1)                      # for reproduceable example
    X <- rdist(1000)                 # sample from X ~ p
    x <- seq(-10,10, .01)
    hist(X, freq=F, breaks=50, xlim=c(-5,5))
    lines(x,p(x),lty=2, col="red")


And when I run the line `dist <-AbscontDistribution(d=p)`, the following error is shown:

```none
Error in seq.default(from = low1, to = upp1, length = ngrid) : 
  'from' cannot be NA, NaN or infinite
Além disso: Warning message:
In sqrt((y - 1)/y + 1) : NaNs produzidos

I'm trying to do something similar to what was done in this link, but with my function. How can I achieve this?

mkrieger1
  • 19,194
  • 5
  • 54
  • 65
david clarck
  • 157
  • 2
  • 11
  • 2
    I plot the graph of your density function, turns out it is not a well defined pdf function: it is not continuous and goes to infinity at multiple points, I wonder it is integrable or not, and it has both positive and negative values. So maybe before asking how to generate random variables from the distribution, make sure your function is really a valid pdf function. – Consistency May 27 '17 at 01:00
  • 1
    Can you elaborate on where the definition for this "density" came from? Is the support supposed to be the whole real line? If so, it does not satisfy the definition of density, as stated above, but it could be valid for certain supports (in that case, it's the limits of your integral that are wrong). – Chris Haug May 27 '17 at 01:54

0 Answers0