Dear R users,
first I take this opportunity to greet all the R community for your
continuous efforts.
I wrote a function to calculate the pdf and cdf of a custom distribution
(mixed gamma model).
The function is the following:
pmixedgamma3 <- function(y, shape1, rate1, shape2, rate2, prev)
{
density.function<-function(x)
{
shape1=shape1
rate1=rate1
shape2=shape2
rate2=rate2
prev=prev
fun<-prev*((rate1^shape1)/gamma(shape1)*(x^(shape1-1))*exp(-rate1*x))
+ (1-prev)*((rate2^shape2)/gamma(shape2)*(x^(shape2-1))*exp(-rate2*x))
return(fun)
}
den<-density.function(y)
p.mixedgamma<-integrate(f=density.function,lower=0,upper=y)
return(list(density=den,cumdensity=p.mixedgamma$value))
}
Why doesn't this calculate cumdensity (p.mixedgamma) in case x is a vector?
For instance try:
pmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5)
As you can see cumdensity is calculated just for the first x value.
Thanks in advance
Best
Mauro
--
Mauro Rossi
Istituto di Ricerca per la Protezione Idrogeologica
Consiglio Nazionale delle Ricerche
Via della Madonna Alta, 126
06128 Perugia
Italia
Tel. +39 075 5014421
Fax +39 075 5014420
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Parte allegato al messaggio
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20120315/3fd18974/attachment.pl>
I believe the difficulty is that the integrate function isn't vectorized: add this line and you should be good pmixedgamma3 <- Vectorize(pmixedgamma3) Read ? Vectorize for details if you only need to vectorize certain arguments. Michael On Thu, Mar 15, 2012 at 6:08 AM, Mauro Rossi <mauro.rossi at irpi.cnr.it> wrote:> Dear R users, > ? ?first I take this opportunity to greet all the R community for your > continuous efforts. > I wrote a function to calculate the pdf and cdf of a custom distribution > (mixed gamma model). > The function is the following: > > pmixedgamma3 <- function(y, shape1, rate1, shape2, rate2, prev) > ? ? ? ?{ > ? ? ? ?density.function<-function(x) > ? ? ? ? ? ? ? ?{ > ? ? ? ? ? ? ? ?shape1=shape1 > ? ? ? ? ? ? ? ?rate1=rate1 > ? ? ? ? ? ? ? ?shape2=shape2 > ? ? ? ? ? ? ? ?rate2=rate2 > ? ? ? ? ? ? ? ?prev=prev > > ?fun<-prev*((rate1^shape1)/gamma(shape1)*(x^(shape1-1))*exp(-rate1*x)) + > (1-prev)*((rate2^shape2)/gamma(shape2)*(x^(shape2-1))*exp(-rate2*x)) > ? ? ? ? ? ? ? ?return(fun) > ? ? ? ? ? ? ? ?} > ? ? ? ?den<-density.function(y) > ? ? ? ?p.mixedgamma<-integrate(f=density.function,lower=0,upper=y) > ? ? ? ?return(list(density=den,cumdensity=p.mixedgamma$value)) > ? ? ? ?} > > Why doesn't this calculate cumdensity (p.mixedgamma) in case x is a vector? > > For instance try: > pmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5) > > As you can see cumdensity is calculated just for the first x value. > > Thanks in advance > Best > Mauro > > -- > Mauro Rossi > Istituto di Ricerca per la Protezione Idrogeologica > Consiglio Nazionale delle Ricerche > Via della Madonna Alta, 126 > 06128 Perugia > Italia > Tel. +39 075 5014421 > Fax +39 075 5014420 > > > ______________________________________________ > R-help at 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. >
On Mar 15, 2012, at 6:08 AM, Mauro Rossi wrote:> Dear R users, > first I take this opportunity to greet all the R community for > your continuous efforts. > I wrote a function to calculate the pdf and cdf of a custom > distribution (mixed gamma model). > The function is the following: > > pmixedgamma3 <- function(y, shape1, rate1, shape2, rate2, prev) > { > density.function<-function(x) > { > shape1=shape1 > rate1=rate1 > shape2=shape2 > rate2=rate2 > prev=prev > fun<-prev*((rate1^shape1)/gamma(shape1)*(x^(shape1-1))*exp(- > rate1*x)) + (1-prev)*((rate2^shape2)/ > gamma(shape2)*(x^(shape2-1))*exp(-rate2*x)) > return(fun) > } > den<-density.function(y) > p.mixedgamma<-integrate(f=density.function,lower=0,upper=y) > return(list(density=den,cumdensity=p.mixedgamma$value)) > } > > Why doesn't this calculate cumdensity (p.mixedgamma) in case x is a > vector?You did not build the function in a manner that would vectorize, but perhaps the convenience vuntion would help here: Vpmixedgamma3 <- Vectorize(pmixedgamma3)> > For instance try: > pmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5)> Vpmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5) [,1] [,2] density 0.12511 0.002908153 cumdensity 0.5420703 0.9950046> > As you can see cumdensity is calculated just for the first x value. > >-- David Winsemius, MD West Hartford, CT