[R] Double Pareto Log Normal Distribution DPLN
David R Forrest
drf at vims.edu
Thu Nov 14 22:29:08 CET 2013
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>
More information about the R-help
mailing list