R help - Mar 2008 - Rejection sampling to draw from distributions

Dear friends,

Please find below the code that I have employed for a rejection sampler to draw
from asymmetric laplace distributions. I was wondering if this code can be
written more efficiently? Are there more efficient ways of drawing random
numbers from asymmetric laplace distributions??

Thanks in advance for your help and have a great weekend.

Regards

Anup

***********************************
ral <- function(n,p,s=3)
{
rout <- matrix(0,n)
for(i in 1:n)
{
  repeat
    {
      root <- rnorm(1)
      ratio <- dal(val=root, p=p)/(s*dnorm(root))
      alpha <- runif(1, min=0, max=1)
      if(alpha<ratio)
        {break}
      }
   rout[i,] <- root
}
return(rout)}


       
---------------------------------

	[[alternative HTML version deleted]]

On Fri, 14 Mar 2008, Anup Nandialath wrote:
> Dear friends,
>
> Please find below the code that I have employed for a rejection sampler 
> to draw from asymmetric laplace distributions. I was wondering if this 
> code can be written more efficiently? Are there more efficient ways of 
> drawing random numbers from asymmetric laplace distributions??
>
Usually.

It really helps to "provide commented, minimal, self-contained, 
reproducible code" as you are asked, but you have not shown what dal() is.

If dal() is **vectorized**, then a vectorized approach like this will 
often work well:

rout <- matrix(0,n)
rejected <- seq(n)
k <- n
while( k > 0 ){
   rout[ rejected ] <- rnorm( k )
   ratio <- dal(val=rout[ rejected ], p=p)/(s*dnorm(rout[ rejected ]))
   alpha <- runif( k , min=0, max=1 )
   rejected <- rejected[ alpha >= ratio ]
   k <- length( rejected )
}

Obviously, I cannot test this without dal().

HTH,

Chuck

p.s. Putting a limit on iterations is usually a good idea, too.

> Thanks in advance for your help and have a great weekend.
>
> Regards
>
> Anup
>
> ***********************************
> ral <- function(n,p,s=3)
> {
> rout <- matrix(0,n)
> for(i in 1:n)
> {
>  repeat
>    {
>      root <- rnorm(1)
>      ratio <- dal(val=root, p=p)/(s*dnorm(root))
>      alpha <- runif(1, min=0, max=1)
>      if(alpha<ratio)
>        {break}
>      }
>   rout[i,] <- root
> }
> return(rout)}
>
>
>
> ---------------------------------
>
> 	[[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.
>
Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901