Re: [R] User-defined random variable

From: Matthias Kohl <Matthias.Kohl_at_uni-bayreuth.de>
Date: Sun 01 May 2005 - 18:45:31 EST

Achim Zeileis wrote:

>On Sat, 30 Apr 2005, Peter Dalgaard wrote:
>
>
>
>>Paul Smith <phhs80@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@biostat.ku.dk) FAX: (+45) 35327907
>>
>>______________________________________________
>>R-help@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@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)

R-help@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 Received on Sun May 01 18:48:13 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:31:31 EST