Search code examples
rdistributionresampling

how to fit a von mises distribution to my data for generating random samples


My data is comprised of 16 pairs of distance and bearings from a particular location. I am trying to generate 1000 resamples of those 16 pairs (i.e. create new sets of X2 and Y2). So that in the end I will have 1000, 16 pairs of distances and bearings resulting in new 16 spatial points.

my data, using the bearing and the distance to generate X2 and Y2

What I have done so far is resample (reshuffle) from the 16 values I already have, but that was no go with my advisors.

f2 <- function(x) data.frame(bearing = sample(min(HRlog$beartoenc):max(HRlog$beartoenc), 16, replace = TRUE), 
                                distance = sample(min(HRlog$distoenc):max(HRlog$distoenc), 16, replace = TRUE))
    
    se1randcent <- as.data.frame(lapply(seq(1000), f2))

I have been told I should resample according to the von mises distribution, i.e fit the distribution to my data then regenerate the 16 pairs from this distribution according the K value I get. I don't really know what this means. Can anyone help me figure it out?


Solution

  • The package circular in R could be helpful. The von Mises distribution kappa parameter can be calculated from the angles you provided either using a minimization method, as I show below, or through built in maximum likelihood estimators mle.vonmises(). Once you have the parameters you can use rvonmises with the number of samples and calculated parameters to generate the sample. The generated sample seems to be on [0,2pi], so there could be some adjustment to make sure the mean values are correctly represented.

    Fitting the distance would probably be a separate distribution and the issue of the possible dependency between the two is not addressed in this answer.

    library(circular) # circular statistics and bessel functions
    
    # converting the bearing to be on the interval [-pi,pi] which is conventional for von Mises
    bearing <- c(19.07,71.88,17.23,202.39,173.67,357.04,5.82,5.82,95.53,5.82,94.13,157.67,19.07,202.39,173.67,128.15)
    bearing_rad <- bearing*2*pi/360 - pi
    
    # sample statistics
    circ_mean <-  mean.circular(bearing_rad) # mu of von Mises
    circ_sd <- sd.circular(bearing_rad) # related to kappa of von Mises
    circ_var <- var.circular(bearing_rad)
    
    # function to return difference in variances between
    diff_vars2 <- function(kappa){
      
      # squaring to make the function convex
      return((1 - A1(kappa) - circ_var)^2)
    }
    
    # solving for kappa by matching the variances
    kappa_solution <- optim(par = 1,fn=diff_vars2,lower = 0,method="L-BFGS-B")
    
    # sample from von mises distribution
    sampled_vals <- rvonmises(n=100, mu=circ_mean, kappa=kappa_solution$par)
    

    Added content based on comments

    One problem with tests for uniformity are that you have a small sample size. Two methods that seem appropriate are the Rayleigh and Kuiper Tests, which test against uniformity. Background on those are given at NCSS Manual

    Both are implemented in circular, but I am not sure if the modified Rayleigh is used. The results for bearings_rad indicate that Rayleigh p-value = 0.2 and Kuiper p-value < 0.05.

    rayleigh.test(x=bearing_rad)
    kuiper.test(x=bearing_rad)
    

    You can add the fitted histogram to the above plot by using dvonmises. This will give the radius, which can be converted to x and y using the standard polar coordinate translation. Making the angles work can be a bit tricky. If you don't want the rose diagram in the background you can use plot.

    rose.diag(bearing_rad)
    
    density_vals <- dvonmises(x=seq(0,2*pi,0.01)-circ_mean,mu = 0,kappa=kappa_solution$par)
    
    x_from_polar <- density_vals*cos(seq(0,2*pi,0.01))
    y_from_polar <- density_vals*sin(seq(0,2*pi,0.01))
    lines(x=x_from_polar,y=y_from_polar,col='red')