Guys 
Any help in this please.
Regards
Bander



> On 14 Nov 2013, at 09:29 pm, "David R Forrest" <drf@vims.edu> wrote:
> 
> 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@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@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@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@vims.edu
>>> To: cs_2002@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@vims.edu > > To: cs_2002@hotmail.com > > CC: r-help@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@vims.edu > > >> To: cs_2002@hotmail.com > > >> CC: r-help@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@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@vims.edu > > >> > > >> > > >> > > >> > > >> > > > > > > [[alternative HTML version deleted]] > > > > > > ______________________________________________ > > > R-help@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@vims.edu > > > > > > > > -- Dr. David Forrest drf@vims.edu 804-684-7900w 757-968-5509h 804-413-7125c 104 Three Point Ct Yorktown, VA, 23692-4325
>> 
>> --
>> Dr. David Forrest
>> drf@vims.edu
> 
> --
> Dr. David Forrest
> drf@vims.edu
> 
> 
> 

	[[alternative HTML version deleted]]

