[R] Double Pareto Log Normal Distribution DPLN

[David R Forrest]

...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" <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)
}

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" <cs_2002 at hotmail.com>
 wrote:

> 
> Hi
> 
> I found this paper http://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" <cs_2002 at hotmail.com>
>> 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




-------------- 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/20131113/37731731/attachment.bin>