[R] User-defined random variable

[Matthias Kohl]

Achim Zeileis wrote:

>On Sat, 30 Apr 2005, Peter Dalgaard wrote:
>
>  
>
>>Paul Smith <phhs80 at gmail.com> writes:
>>
>>    
>>
>>>Dear All
>>>
>>>I would like to know whether it is possible with R to define a
>>>discrete random variable different from the ones already defined
>>>inside R and generate random numbers from that user-defined
>>>distribution.
>>>      
>>>
>>Yes. One generic way is to specify the quantile function (as in
>>qpois() etc.) and do qfun(runif(N)).
>>    
>>
>
>If the support discrete but also finite, you can also use sample(), e.g.
>
>  sample(myset, N, replace = TRUE, prob = myprob)
>
>Z
>
>  
>
>>--
>>   O__  ---- Peter Dalgaard             Blegdamsvej 3
>>  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N
>> (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
>>~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
>>
>>______________________________________________
>>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
>>
>>    
>>
>
>______________________________________________
>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
>  
>
Hi,

one can also use our R package "distr" to generate discrete random 
variables. The subsequent code provides a function which generates an 
object of class "DiscreteDistribution" based on a finite support "supp". 
If "prob" is missing all elements in "supp" are equally weighted.

Matthias


# generating function
DiscreteDistribution <- function(supp, prob){
    if(!is.numeric(supp))
        stop("'supp' is no numeric vector")
    if(any(!is.finite(supp)))   # admit +/- Inf?
        stop("inifinite or missing values in supp")
    len <- length(supp)
    if(missing(prob)){
        prob <- rep(1/len, len)
    }else{
        if(len != length(prob))
            stop("'supp' and 'prob' must have equal lengths")
        if(any(!is.finite(prob)))
            stop("inifinite or missing values in prob")
        if(!identical(all.equal(sum(prob), 1,
                            tolerance = 2*distr::TruncQuantile), TRUE))
            stop("sum of 'prob' has to be (approximately) 1")
        if(!all(prob >= 0))
            stop("'prob' contains values < 0")
    }
    if(length(usupp <- unique(supp)) < len){
        warning("collapsing to unique support values")
        prob <- as.vector(tapply(prob, supp, sum))
        supp <- sort(usupp)
        len <- length(supp)
    }else{
        o <- order(supp)
        supp <- supp[o]
        prob <- prob[o]
    }
   
    if(len > 1){
      if(min(supp[2:len] - supp[1:(len - 1)]) < distr::DistrResolution)
        stop("grid too narrow --> change DistrResolution")
    }
    rfun <- function(n){
        sample(x = supp, size = n, replace = TRUE, prob = prob)
    }
    intervall <- distr::DistrResolution / 2 
 
    supp.grid <- as.numeric(matrix(
                      rbind(supp - intervall, supp + intervall), nrow = 1))
    prob.grid <- c(as.numeric(matrix(rbind(0, prob), nrow = 1)), 0)
    dfun <- function(x){ stepfun(x = supp.grid, y = prob.grid)(x) }
 
    cumprob <- cumsum(prob)
    pfun <- function(x){ stepfun(x = supp, y = c(0, cumprob))(x) }

    qfun <- function(x){ supp[sum(cumprob<x)+1] }

    return(new("DiscreteDistribution", r = rfun, d = dfun, p = pfun,
        q = qfun, support = supp))
}

# some examples
supp <- rpois(20, lambda=7) # some support vector
D1 <- DiscreteDistribution(supp = supp) # prob is missing
r(D1)(10) # 10 random numbers of Distribution D1
support(D1) # support
d(D1)(support(D1)) # pdf
p(D1)(5) # cdf
q(D1)(0.5) # quantile (here: median)
plot(D1) # plot of pdf, cdf and quantile
D2 <- lgamma(D1) # apply member of group generic "Math"
plot(D2)
D3 <- D1 + Binom(size=10) # convolution with object of class 
"DiscreteDistribution"
plot(D3)
D4 <- D1 + Norm() # convolution with object of class "AbscontDistribution"
plot(D4)