[R] Plot cumulative probability of beta-prime distribution

David Winsemius dwinsemius at comcast.net
Thu Jul 2 14:43:07 CEST 2009


I am not at all familiar with that distribution but the obvious  
solution would appear to be:

?integrate

 > BetaprimeDensity <- function(x) x^(shape1-1) * (1+x)^(-shape1- 
shape2) / beta(shape1,shape2)
 > shape1 <- 1
 > shape2 <-1
 > integrate(BetaprimeDensity, 0 , 1)
0.5 with absolute error < 5.6e-15

If I knew more about that distribution I would be able to see if that  
were sensible, but severe caveats apply, since I just took you  
expression and massaged it into a function. If I knew what the domain  
was supposed to be, it would be simple matter to wrap that expression  
into a function that could be called pbetaprime. Perhaps:

 > pbetaprime <- function(ul, ll=0,  s1=2, s2=2) integrate(function(x)  
{x^(s1-1) * (1+x)^(-s1-s2) / beta(s1,s2)}, ll, ul)
 > pbetaprime(5, , 5, 3)
0.9042245 with absolute error < 7e-05

Does that "work" when given the proper sequence vector of values and  
submitted to plot?

-- 
DW


On Jul 1, 2009, at 3:19 AM, aledanda wrote:

>
> Hallo,
>
> I need your help.
> I fitted my distribution of data with beta-prime, I need now to plot  
> the
> Cumulative distribution. For other distribution like Gamma is easy:
>
> x <- seq (0, 100, 0.5)
> plot(x,pgamma(x, shape, scale), type= "l", col="red")
>
> but what about beta-prime? In R it exists only pbeta which is  
> intended only
> for the beta distribution (not beta-prime)
>
> This is what I used for the estimation of the parameters:
>
> mleBetaprime <- function(x,start=c(1,1)) {
>   mle.Estim <- function(par) {
>     shape1 <- par[1]
>     shape2 <- par[2]
>     BetaprimeDensity <- NULL
>     for(i in 1:length(posT))
>       BetaprimeDensity[i] <- posT[i]^(shape1-1) *
> (1+posT[i])^(-shape1-shape2) / beta(shape1,shape2)
>     return(-sum(log(BetaprimeDensity)))
>    }
>   est <- optim(fn=mle.Estim,par=start,method="Nelder-Mead")
>   return(list(shape1=est$par[1],shape2=est$par[2]))
>  }
> posbeta1par <- fdp(posT, family= "beta1")
>
> Hope you can help me.
>
> Thanks a lot!!!

David Winsemius, MD
Heritage Laboratories
West Hartford, CT




More information about the R-help mailing list