[R] simulating from Langevin distributions

[Ranjan Maitra]

Thanks to Ravi Varadhan for providing the solution. I guess the answer is that if 

X is MVN with mean mu and dispersion matrix given by I/K, then X/norm(X) is Langevin with the required parameters. A reference for this is Watson's Statistics on Spheres.

Many thanks again and best wishes.
Ranjan


On Tue, 13 Feb 2007 14:44:08 -0500 "Ravi Varadhan" <rvaradhan at jhmi.edu> wrote:

> Hi Ranjan,
> 
> I think that the following would work:
> 
> library(MASS)
> 
> rlangevin <- function(n, mu, K) {
> q <- length(mu)
> norm.sim <- mvrnorm(n, mu=mu, Sigma=diag(1/K, q))
> cp <- apply(norm.sim, 1, function(x) sqrt(crossprod(x)))
> sweep(norm.sim, 1, cp, FUN="/")
> }
> 
> > mu <- runif(7)
> > mu <- mu / sqrt(crossprod(mu))
> > K <- 1.2
> > ylang <- rlangevin(n=10, mu=mu, K=K)
> > apply(ylang,1,crossprod)
>  [1] 1 1 1 1 1 1 1 1 1 1
> >
> 
> I hope that this helps.
> 
> Ravi.
> ----------------------------------------------------------------------------
> -------
> 
> Ravi Varadhan, Ph.D.
> 
> Assistant Professor, The Center on Aging and Health
> 
> Division of Geriatric Medicine and Gerontology 
> 
> Johns Hopkins University
> 
> Ph: (410) 502-2619
> 
> Fax: (410) 614-9625
> 
> Email: rvaradhan at jhmi.edu
> 
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
>  
> 
> ----------------------------------------------------------------------------
> --------
> 
> 
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Ranjan Maitra
> Sent: Tuesday, February 13, 2007 1:04 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] simulating from Langevin distributions
> 
> Dear all,
> 
> I have been looking for a while for ways to simulate from Langevin
> distributions and I thought I would ask here. I am ok with finding an
> algorithmic reference, though of course, a R package would be stupendous!
> 
> Btw, just to clarify, the Langevin distribution with (mu, K), where mu is a
> vector and K>0 the concentration parameter is defined to be:
> 
> f(x) = exp(K*mu'x) / const where both mu and x are p-variate vectors with
> norm 1.
> 
> For p=2, this corresponds to von-Mises (for which algorithms exist,
> including in R/Splus) while for p=3, I believe it is called the Fisher
> distribution. I am looking for general p.
> 
> Can anyone please help in this?
> 
> Many thanks and best wishes,
> Ranjan
> 
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>