...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 protected]> 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" <[email protected]>
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: [email protected]
>> To: [email protected]
>> CC: [email protected]
>> 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" <[email protected]>
>> 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]]
>>>
>>> ______________________________________________
>>> [email protected] 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
>> [email protected]
>>
>>
>>
>>
>>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [email protected] 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
[email protected]
signature.asc
Description: Message signed with OpenPGP using GPGMail
______________________________________________ [email protected] 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.
