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 >> 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: d...@vims.edu > > >> To: cs_2...@hotmail.com > > >> CC: >> r-h...@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 > > >> 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
signature.asc
Description: Message signed with OpenPGP using GPGMail
______________________________________________ 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.
