conjugate
The simplest way to use the Von-Mises distribution in pcvr is through the conjugate function, where “vonmises” and “vonmises2” are valid methods. As with other conjugate methods these are implemented for single or multi value traits, but unlike other methods these are only necessarily supported in comparing to other samples from the same distribution. If more distributions become tenable to add as circular or wrapped functions then this may be revisted.
“vonmises” method
The “vonmises” method uses the fact that the conjugate prior for the direction parameter (\(\mu\)) is itself a Von-Mises distribution. Utilizing this conjugacy requires that we assume a known \(\kappa\) for the complete distribution so that updating the \(\mu\) parameter is straightforward. Conceptually it may be helpful to consider this similarly to the “T” method for comparing the means of guassians.
Priors for this method should specify a list containing “mu”, “kappa”, “boundary”, “known_kappa”, and “n” elements. In that prior “mu” is the direction of the circular distribution, “kappa” is the precision of the mean, “boundary” is a vector including the two values that are the where the circular data “wraps” around the circle, “known_kappa” is the fixed value of precision for the total distribution, and “n” is the number of prior observations. If the prior is not specified then the default is list(mu = 0, kappa = 1, boundary = c(-pi, pi), known_kappa = 1, n = 1). As per other methods for the conjugate function, the “posterior” part of the output is of the same form as the prior.
Example
First we’ll simulate some multi value data
Code
mv_gauss <- mvSim(
dists = list(
rnorm = list(mean = 50, sd = 10),
rnorm = list(mean = 60, sd = 12)
),
n_samples = c(30, 40)
)Next we’ll run conjugate specifying that our data is on a circle defined over [0, 180] with an expected direction around 45 (90 degrees on the full [0,360] or \(\pi/2\) radians) and low precision.
Code
vm_ex1 <- conjugate(
s1 = mv_gauss[1:30, -1],
s2 = mv_gauss[31:70, -1],
method = "vonmises",
priors = list(mu = 45, kappa = 1, boundary = c(0, 180), known_kappa = 1, n = 1),
plot = TRUE, rope_range = c(-5, 5), rope_ci = 0.89,
cred.int.level = 0.89, hypothesis = "equal"
)The summary shows normal conjugate output, here showing a posterior probability of ~91% that our samples have equal means (remember the difference in our simulated data is now on a circle).
Code
vm_ex1$summary## HDE_1 HDI_1_low HDI_1_high HDE_2 HDI_2_low HDI_2_high hyp post.prob
## 1 49.46678 32.66822 64.46699 59.52812 41.66264 79.59501 equal 0.6375772
## HDE_rope HDI_rope_low HDI_rope_high rope_prob
## 1 -10.49173 -34.87026 14.28739 0.2237951Displaying plots of these data can be slower than for other conjugate methods due to the density of the support. To explain, the Von-Mises distribution is defined in on the unit circle [\(-\pi\), \(\pi\)] so in order to have support that works to project that data into whatever space the boundary in the prior specifies the support has to be very dense.
Note also that our rope_range is specified in the boundary units space, which is not necessarily the unit circle.
Code
vm_ex1$plot # not printed due to being a very dense ggplotWe get very similar results using roughly analogous single value traits.
Code
vm_ex1_1 <- conjugate(
s1 = rnorm(30, 50, 10),
s2 = rnorm(40, 60, 12),
method = "vonmises",
priors = list(mu = 0, kappa = 1, known_kappa = 1, boundary = c(0, 180), n = 1),
plot = FALSE, rope_range = c(-0.1, 0.1), rope_ci = 0.89,
cred.int.level = 0.89, hypothesis = "equal"
)
do.call(rbind, vm_ex1_1$posterior)## mu kappa known_kappa n boundary
## [1,] 50.04425 7.002869 1 31 numeric,2
## [2,] 57.94215 8.385856 1 41 numeric,2Single value traits work in the same way. Note that if we omit parts of the prior then they will be filled in with the default prior values.
Code
set.seed(42)
vm_ex2 <- conjugate(
s1 = brms::rvon_mises(100, -3.1, 2),
s2 = brms::rvon_mises(100, 3.1, 2),
method = "vonmises",
priors = list(mu = 0, kappa = 1, known_kappa = 2),
plot = TRUE, rope_range = c(-0.1, 0.1), rope_ci = 0.89,
cred.int.level = 0.89, hypothesis = "equal"
)We check our summary and see around 75% chance that these are equal
Code
vm_ex2$summary## HDE_1 HDI_1_low HDI_1_high HDE_2 HDI_2_low HDI_2_high hyp post.prob
## 1 -3.071143 2.393966 -2.224163 3.140433 2.381274 -2.297627 equal 0.940651
## HDE_rope HDI_rope_low HDI_rope_high rope_prob
## 1 0.07243397 -1.081833 1.232842 0.1241434Code
do.call(rbind, vm_ex2$posterior)## mu kappa known_kappa n boundary
## [1,] -3.071143 4.204577 2 101 numeric,2
## [2,] 3.140433 4.499838 2 101 numeric,2Here our plot is much faster to make since the support is a roughly a thirtieth the size of the previous example.
Code
vm_ex2$plot # not printed due to being a very dense ggplotSometimes it may be helpful to use polar coordinates to consider this data, although limitations in plotting area style geometries in the polar coordinates can be frustrating.
Code
p <- vm_ex2$plot
p[[1]] <- p[[1]] +
ggplot2::coord_polar() +
ggplot2::scale_y_continuous(limits = c(-pi, pi))“vonmises2” method
The “vonmises2” method updates \(\mu\) and \(\kappa\) of the complete Von-Mises distribution. This is done by first taking a weighted average of the prior \(\kappa\) and the MLE of \(\kappa\) based on the sample data then updating \(\mu\) as above.
Priors for this method should specify “mu”, “kappa”, “boundary”, and “n”. Where “mu” is still mean direction, “kappa” is the precision, and boundary/n are as above.
Example
Using the same test data as above we can run the “vonmises2” method.
Code
vm2_ex1 <- conjugate(
s1 = mv_gauss[1:30, -1],
s2 = mv_gauss[31:70, -1],
method = "vonmises2",
priors = list(mu = 45, kappa = 1, boundary = c(0, 180), n = 1),
plot = TRUE, rope_range = c(-5, 5), rope_ci = 0.89,
cred.int.level = 0.89, hypothesis = "equal"
)Code
do.call(rbind, vm2_ex1$posterior)## mu kappa n boundary
## [1,] 49.46678 8.371329 31 numeric,2
## [2,] 59.52812 6.07 41 numeric,2Code
set.seed(42)
vm2_ex2 <- conjugate(
s1 = brms::rvon_mises(100, -3.1, 2),
s2 = brms::rvon_mises(100, 3.1, 2),
method = "vonmises2",
priors = list(mu = 0, kappa = 1),
plot = TRUE, rope_range = c(-0.75, 0.75), rope_ci = 0.89,
cred.int.level = 0.89, hypothesis = "equal"
)Code
vm2_ex2$plot # much lighter to print this since it is in radians
Code
do.call(rbind, vm2_ex2$posterior)## mu kappa n boundary
## [1,] -3.071143 2.091375 101 numeric,2
## [2,] 3.140433 2.237544 101 numeric,2