0

I wish to apply a custom function to each element of a matrix whilst also using elements of a different matrix as inputs to the function.

Specifically, my function generates random samples from a von Mises distribution (circular normal distribution), calling the Rfast package's rvonmises function.

I have one matrix (radians) which records the angle I wish to use for the central tendency of the random generation (similar to the mean), and another matrix (kappa) which records the concentration parameter of the von Mises I wish to use (similar to standard deviation).

I wish to use (for example) element [1, 1] of the radians matrix together with element [1, 1] of the kappa matrix in a call to the von Mises random generator. So, my call for one element would be:

rvonmises(n = 1, m = radians[1, 1], k = kappa[1, 1])

But of course I want this applied across all elements of the matrices. (The rvonmises function doesn't accept multiple m or k values, so for example I couldn't use rvonmises(4, m = c(1, 2, 3, 4), k = c(1, 1.2, 1.4, 1.6)).)

To summarise: I am basically after a more principled (and faster!) way of doing this:

for(i in 1:nrow(radians)){
  for(j in 1:ncol(radians)){
    result[i, j] <- Rfast::rvonmises(1, radians[i, j], kappa[i, j])
  }
}

What I have tried

Based on this post, I have tried to use mapply:

library(Rfast)
set.seed(42)

# random radians to use as input
radians <- matrix(data = runif(12, 0, 2 * pi), 
                  ncol = 4)

# random concentration parameters of the von Mises distribution
kappa <- matrix(data = rgamma(12, 70, 30), 
                ncol = 4)

# function to generate random von mises sample with angle x and 
# concentration parameter k
my_function <- function(m, k){
  Rfast::rvonmises(1, m, k)
}

# my attempt
out <- matrix(mapply(my_function, m = as.data.frame(radians), k = kappa), 
              ncol = 4, byrow = TRUE)

However, I don't think this is working. For example, if I test it by the following (where the central tendency in test_radians increases steadily and I use large values for kappa which leads to precise estimates):

test_radians <- matrix(data = seq(from = 1, to = 2 * pi, length.out = 12), 
                       ncol = 4)

test_kappa <- matrix(data = rep(20, times = 12), 
                     ncol = 4)

test <- matrix(mapply(my_function, m = as.data.frame(test_radians), 
                      k = test_kappa), 
               ncol = 4, byrow = TRUE)

test[1, 1] should be smaller (on average), and test[3, 4] should be largest. (I know due to random variability this won't always be the case, but I've tried it with many replications.)

So, the mapping and matching between matrices isn't working as I had anticipated.

Any guidance welcomed.

JimGrange
  • 262
  • 2
  • 10
  • I'm not familiar with `Rfast::rvonmises`, but running `Rfast::rvonmises(100, 0, 1)` returns this histogram: https://i.imgur.com/p8FurOa.png showing that many values can end up far beyond 0. – Phil Mar 14 '21 at 16:47
  • Yes, the distribution is bound between 0 and 2*pi (see the data in my **radians** matrix). – JimGrange Mar 14 '21 at 16:49
  • RIght, but you say `test[1, 1]` should be around zero, but that's not my experience. – Phil Mar 14 '21 at 16:50
  • Sorry, I now see what you mean. I guess I just need to confirm whether the mapply call is doing what I think it is, but I can't seem to verify it. – JimGrange Mar 14 '21 at 16:55
  • OK I've edited my "test" part a little. I've now used a large value for kappa in the RNG. This produces stable estimates, but this is not what is seen in the resulting test_out matrix. (For example, ```mean(rvonmises(10000, 4, 20))``` gives a result very close to expected mean value (4). This is not seen in the test_out matrix suggesting the mapping is not working as expected.) – JimGrange Mar 14 '21 at 17:12
  • If your concern is just really on whether you are using `mapply()` properly, I would just change your custom function to something simple like `function(m, k) { m * k }`, and then check that the values across the matrix should be what you expect. – Phil Mar 14 '21 at 17:47
  • Thanks, but my question is here because it is **not** working, and I do not understand what error I am making. ```function(m, k) {m * k}``` returns a 9 * 4 matrix, so the ```mapply()``` is not working as I expect it to. – JimGrange Mar 14 '21 at 18:37
  • `matrix(mapply(my_function, m = test_radians, k = test_kappa), ncol = 4, byrow = TRUE)` returns a matrix of 4 columns and 3 rows for me. – Phil Mar 14 '21 at 21:59

1 Answers1

1

You cannot compute the mean of circular observations by simply calling "mean". This is wrong. The correct way is to compute the mean of the cosinus and sinus of the angles and then use the arc tangent. See pcakcges for directional or circular data for this.

Secondly, you gave us an idea, to return a matrix of von Mises generated data. But, since brms does this job for you, at the moment I would go there.

Michail
  • 177
  • 1
  • 6