Hi Bander,
I'm pushing this discussion back to the list, because I'm not sure of
the shape/rate parameters for rpareto and rexp and how they'd be applied
across this mix of typo'd papers.
# Reed Equation 6 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf exponentiated
per end of sec 3:
rdpln<-function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}
# Reed Equation 10:
library(VGAM)
rdpln2<-function(n,a,b,v,t){ rlnorm(n,meanlog=v,sdlog=t)*
rpareto(n,location=1,shape=1/a)/
rpareto(n,location=1,shape=1/b)}
boxplot(data.frame(x1= log(rdpln(a=2.5,b=.01,t=0.45,v=6.5,n=100000)), x2=
log(rdpln2(a=2.5,b=.01,t=0.45,v=6.5,n=100000))))
# Reed equation 8 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf
# with S1 errata #1 from
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964
ddpln <- function(x, a=1, b=1, v=0, t=1){
# Density of Double Pareto LogNormal distribution
# from "b. alzahrani" <cs_2002 at hotmail.com> email of
2013-11-13
# per formula 8 from http://cs.stanford.edu/people/jure/pubs/dpln-kdd08.pdf
c <- (a * b /(a+b))
norm1<-pnorm((log(x)-v-(a*t^2))/t)
norm2<-pnorm((log(x)-v+(b*t^2))/t)
expo1<- a*v+(a^2*t^2/2)
expo2<- -b*v+(b^2*t^2/2) # edited from the paper's eqn 8 with: s/t/v/
z<- (x^(-a-1)) * exp(expo1)*( norm1)
y<- (x^( b-1)) * exp(expo2)*(1-norm2) # 1-norm is the complementary CDF of
N(0,1)
c*(z+y)
}
Dave
On Nov 14, 2013, at 9:12 AM, David Forrest <drf at vims.edu>
wrote:
>
> I think exponentiation of eqn 6 from the Reed paper generates DPLN variates
directly, so maybe:
>
>
rdpln=function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}
>
>
> Dave
>
>
> On Nov 13, 2013, at 4:34 PM, "b. alzahrani" <cs_2002 at
hotmail.com>
> wrote:
>
>> You help is much appreciated. Just one last point if you could help,
Now I want to pass this curve to a function that can generate random numbers
distributed according to DPLN ( right curve).
>>
>> I found the package Runuran can do that
http://cran.r-project.org/web/packages/Runuran/Runuran.pdf but I do not know
how to do it I think it would be something similar to page 8 and 9.
>>
>> Regards
>> ******************************************************************
>> Bander
>> *************************************
>>
>>
>>
>> From: drf at vims.edu
>> To: cs_2002 at hotmail.com
>> Subject: Re: [R] Double Pareto Log Normal Distribution DPLN
>> Date: Wed, 13 Nov 2013 21:13:43 +0000
>>
>>
>> I read the parameters in Fig 4, right as
"DPLN[2.5,0.1,0.45,6.5]", so:
>>
>> x<- 10^seq(0,4,by=.1)
>>
plot(x,ddpln(x,a=2.5,b=.01,v=6.5,t=0.45),log='xy',type='l')
>>
>> ... and the attached graph does not look dissimilar the figure--It
starts at 10^-2, goes through 10^-4 at x=100, and elbows down around 900 and
passes through 10^-6 at about 2000.
>>
>> The correction of Reed helps -- The uncorrected Reed Eq9 equation
suggests that the the 't' in Sehshadri Eq9 should be a 'v' , but
it doesn't exactly make sense with the extra 'a' in there. If the
errata clears that up, then your expo2 term looks just like the expo1 term, but
with a=-b.
>>
>>
>>
>>
>>
>> On Nov 13, 2013, at 3:43 PM, "b. alzahrani" wrote: > Thank
you very much for the help and the change you suggested in my code, I also found
a correction on equation 9 that has been published by Reed ( here
http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf i.e. the original paper on
Double Pareto Log Normal Distribution ). > > can you please see the
correction in this
linkhttp://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964
(in Supporting Information section, appendix S1), does your suggested code
coincide with the correction on this link? as I can see > > Actually, I am
interested in the most right curve in figure 4. and if we plot the curve with
same order of the paper's parameters you suggests:
plot(x,ddpln(x,a=2.8,b=.01,v=6.5,t=0.45),log='xy',type='l') the
curve is different from the one in the paper?? > > Thanks > >
****************************************************************** > Bander
Alzahrani, > ************************************* > > > > >
From: drf at vims.edu > > To: cs_2002 at hotmail.com > > CC: r-help
at stat.math.ethz.ch > > Subject: Re: [R] Double Pareto Log Normal
Distribution DPLN > > Date: Wed, 13 Nov 2013 19:09:34 +0000 > > >
> ...Additionally...your set of parameters match none of the curves in figure
4. > > > > I think the ordering of the parameters as listed on the
graphs is different than in the text of the article. > > > > The
'v' parameter controls the location of the 'elbow' and should be
near log(x) in each graph, while the 't' parameter is the sharpness of
the elbow. So eyeballing the elbows: > > > > log(c(80,150,300))=
4.382027 5.010635 5.703782 > > > > These appear to match positional
parameter #4 in the legends, which you apply as parameter t in your function.
> > > > > > # editing your function a bit: > > > >
ddpln <- function(x, a=2.5, b=0.01, v=log(100), t=v/10){ > > # Density
of Double Pareto LogNormal distribution > > # from "b.
alzahrani" email of 2013-11-13 > > # per formula 8 from
http://cs.stanford.edu/people/jure/pubs/dpln-kdd08.pdf > > > > c
<- (a * b /(a+b)) > > > > norm1<-pnorm((log(x)-v-(a*t^2))/t)
> > norm2<-pnorm((log(x)-v+(b*t^2))/t) > > expo1<-
a*v+(a^2*t^2/2) > > expo2<- -b*v+(b^2*t^2/2) # edited from the
paper's eqn 8 with: s/t/v/ > > > > z<- (x^(-a-1)) *
exp(expo1)*(norm1) > > y<- (x^(b-1)) * exp(expo2)*(1-norm2) # 1-norm is
the complementary CDF of N(0,1) > > > > c*(z+y) > > } >
> > > x<-10^seq(0,5,by=0.1) # 0-10k > > > >
plot(x,ddpln(x,a=2.8,b=.01,v=3.8,t=0.35),log='xy',type='l') #
Similar to fig 4 left. > > > >
plot(x,ddpln(x,a=2.5,b=.01,v=log(2)),log='xy',type='l') >
> plot(x,ddpln(x,a=2.5,b=.01,v=log(10)),log='xy',type='l')
> >
plot(x,ddpln(x,a=2.5,b=.01,v=log(100)),log='xy',type='l') >
> plot(x,ddpln(x,a=2.5,b=.01,v=log(1000)),log='xy',type='l')
> >
plot(x,ddpln(x,a=2.5,b=.01,v=log(10000)),log='xy',type='l') >
> > > Dave > > > > On Nov 13, 2013, at 11:43 AM, "b.
alzahrani" > > wrote: > > > > > > > > Hi >
> > > > > I found this
paperhttp://cs.stanford.edu/people/jure/pubs/dpln-kdd08.pdf that models the DPLN
distribution as in equation 8. I implemented this in R but cannot get the same
curve as in Figure 4. can you please check if my code below is correct: e.g. is
the use of pnorm() correct here? > > > > > > ddlpn <-
function(x){ > > > a=2.8 > > > b=0.01 > > > v=0.45
> > > t=6.5 > > > j <- (a * b /(a+b)) > > > >
> > norm1<-pnorm((log(x)-v-(a*t^2))/t) > > > expo1<-
a*v+(a^2*t^2/2) > > > > > >
z<-exp(expo1)*(x^(-a-1))*(norm1) > > > > > >
norm2<-pnorm((log(x)-v+(b*t^2))/t) > > > expo2<- -b*t+(b^2*t^2/2)
> > > > > > y<- x^(b-1)*exp(expo2)*(1-norm2) # 1-norm is
the complementary CDF of N(0,1) > > > j*(z+y) > > > } >
> > ******************************************************************
> > > Bander Alzahrani, Teacher Assistant > > > Information
Systems Department > > > Faculty of Computing & Information
Technology > > > King Abdulaziz University > > > > >
> ************************************* > > > > > > >
> > > > >> From: drf at vims.edu > > >> To:
cs_2002 at hotmail.com > > >> CC: r-help at stat.math.ethz.ch >
> >> Subject: Re: [R] Double Pareto Log Normal Distribution > >
>> Date: Tue, 12 Nov 2013 16:51:22 +0000 > > >> > >
>> > > >> http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf is
the original ref and has the equations. > > >> > > >>
library(VGAM) for *pareto() and library(stats) for *lnorm() should get you most
of the way there. > > >> > > >> On Nov 12, 2013, at
10:47 AM, "b. alzahrani" > > >> wrote: > > >>
> > >>> Hi guys > > >>> I would like to generate
random number Double Pareto Log Normal Distribution (DPLN). does anyone know how
to do this in R or if there is any built-in function. > > >>>
> > >>> Thanks > > >>> > > >>>
****************************************************************** > >
>>> Bander > > >>> *************************************
> > >>> > > >>> > > >>> > >
>>> [[alternative HTML version deleted]] > > >>> >
> >>> ______________________________________________ > >
>>> 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. > > >> > > >> -- > > >>
Dr. David Forrest > > >> drf at vims.edu > > >> >
> >> > > >> > > >> > > >> > >
> > > > [[alternative HTML version deleted]] > > > >
> > ______________________________________________ > > > 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. > > >
> -- > > Dr. David Forrest > >drf at vims.edu > > > >
> > > > -- Dr. David Forrest drf at vims.edu 804-684-7900w
757-968-5509h 804-413-7125c 104 Three Point Ct Yorktown, VA, 23692-4325
>
> --
> Dr. David Forrest
> drf at vims.edu
>
--
Dr. David Forrest
drf at vims.edu
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20131114/354f988a/attachment.bin>