[R] simulating from Langevin distributions

[Ravi Varadhan]

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.