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

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to