Thanks Dave for your help on this. Waiting other suggestions from the list on this.
Regards Bander > On 14 Nov 2013, at 09:29 pm, "David R Forrest" <d...@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_2...@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 <d...@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_2...@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: d...@vims.edu >>> To: cs_2...@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: d...@vims.edu > > To: cs_2...@hotmail.com > > CC: > > r-h...@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 ht! tp://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 D! PLN 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: d...@vims.edu > > >> To: cs_2...@hotmail.com > > >> CC:! r-h...@stat.math.ethz.ch > > >> Subject: Re: [R] Double Pareto Log No rmal 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 > > >> d...@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 > >d...@vims.edu > > > > > > > > -- Dr. David Forrest d...@vims.edu 804-684-7900w 757-968-5509h 804-413-7125c 104 Three Point Ct Yorktown, VA, 23692-4325 >> >> -- >> Dr. David Forrest >> d...@vims.edu > > -- > Dr. David Forrest > d...@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.
