# [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.
>>
>>
>>
>>
>>
>
> --
> 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>
```