| Title: | Circular Statistics, from "Topics in Circular Statistics" (2001) |
|---|---|
| Description: | Circular Statistics, from "Topics in Circular Statistics" (2001) S. Rao Jammalamadaka and A. SenGupta, World Scientific. |
| Authors: | S-plus original by Ulric Lund <[email protected]>, R port by Claudio Agostinelli <[email protected]> |
| Maintainer: | Claudio Agostinelli <[email protected]> |
| License: | GPL-2 |
| Version: | 0.2-6 |
| Built: | 2024-11-20 06:24:04 UTC |
| Source: | CRAN |
Evaluates the first and zeroth order Bessel functions of the first kind at a specified non-negative real number, and returns the ratio.
A1(kappa)A1(kappa)
kappa |
non-negative numeric value at which to evaluate the Bessel functions. |
The function use besselI.
If I1(kappa) is the first order Bessel function and I0(kappa) is the zeroth order Bessel function, then A1(kappa) returns I1(kappa)/I0(kappa).
Inverse function of the ratio of the first and zeroth order Bessel functions of the first kind. This function is used to compute the maximum likelihood estimate of the concentration parameter of a von Mises distribution.
A1inv(x)A1inv(x)
x |
numeric value in the interval between 0 and 1. |
A1inv(0) = 0 and A1inv(1) = inf. This function is useful in estimating the concentration parameter of data from a von Mises distribution.
Returns the value k, such that A1inv(x) = k, i.e. A1(k) = x.
#Generate data from a von Mises distribution data <- rvm(50, pi, 4) #Estimate the concentration parameter s <- sum(sin(data)) c <- sum(cos(data)) mean.dir <- atan2(s, c) kappa <- A1inv(mean(cos(data - mean.dir)))#Generate data from a von Mises distribution data <- rvm(50, pi, 4) #Estimate the concentration parameter s <- sum(sin(data)) c <- sum(cos(data)) mean.dir <- atan2(s, c) kappa <- A1inv(mean(cos(data - mean.dir)))
Tests for a change in mean direction, concentration, or both, given a set of directional data points.
change.pt(x)change.pt(x)
x |
vector of angular measurements in radians. |
In either context, the user can choose which statistic (max or ave) to use, and then consult the appropriate table provided in the book referenced below. The critical values for these 4 statistics are to be found in Table 11.3 (or Figure 11.3) for rmax, Table 11.4 (or Figure 11.4) for rave, Figure 11.5 for tmax and Figure 11.6 for tave.
Returns a data frame with variables n, rho=r/n, rmax, k.r, rave, tmax, k.t, and tave. The first of these is the sample size, followed by the overall mean resultant length. Both of these are needed to enter any of the tables or nomograms(see under Details). The other values represent the change point test statistics. While rmax and rave test for a change in mean direction (with unknown concentration), tmax and tave are useful in the context of testing more generally, for a change in mean direction and/or concentration. k.r and k.t are the observation numbers for which rmax and tmax attain their maximum value and indicate the observation at which the change is most likely to have occurred, when the tables or nomograms indicate significance.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Chapter 11, World Scientific Press, Singapore.
Computes a circular version of the Pearson's product moment correlation, and performs a significance test if requested.
circ.cor(alpha, beta, test=FALSE)circ.cor(alpha, beta, test=FALSE)
alpha |
vector of circular data measured in radians. |
beta |
vector of circular data measured in radians. |
test |
if test = TRUE, then a significance test for the correlation coefficient is computed. |
The correlation coefficient is computed like Pearson's product moment correlation for two linear variables X and Y. In the computational formula, however, (xi - xbar) and (yi - ybar) are replaced by sin(xi - xbar) and sin(yi - ybar), where xbar and ybar in the second two expressions are the mean directions of the samples.
Returns a data frame with variables r, a circular version of the Pearson's product moment correlation, test.stat and p.value, the test statistic and p-value respectively, for testing significance of the correlation coefficient. test.stat and p.value are by default not produced, but are given when test=TRUE is specified in the function call.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 8.2, World Scientific Press, Singapore.
Jammalamadaka, S. and Sarma, Y. (1988). A correlation coefficient for angular variables. Statistical Theory and Data Analysis 2. North Holland: New York.
# Generate two circular data sets, and compute their correlation. data1 <- rvm(50, 0, 3) data2 <- data1 + pi + rvm(50, 0, 10) circ.cor(data1, data2, test=TRUE)# Generate two circular data sets, and compute their correlation. data1 <- rvm(50, 0, 3) data2 <- data1 + pi + rvm(50, 0, 10) circ.cor(data1, data2, test=TRUE)
Computes measures of dispersion for a directional data set.
circ.disp(x)circ.disp(x)
x |
vector of circular data measured in radians. |
Returns a dataframe with the following components. The sample size, n; the resultant length, r; the mean resultant length, rbar= r/n; and the circular variance, var=(1-r/n).
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 1.3, World Scientific Press, Singapore.
circ.mean, circ.summary, est.kappa, est.rho.
Returns the mean direction of a vector of circular data.
circ.mean(x)circ.mean(x)
x |
vector of data points measured in radians. |
Each observation is treated as a unit vector, or point on the unit circle. The resultant vector of the observations is found, and the direction of the resultant vector is returned.
Returns the mean direction of the data.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 1.3, World Scientific Press, Singapore.
circ.disp, circ.summary, est.kappa, est.rho.
# Compute the mean direction of a random sample of observations. data <- runif(50, 0, pi) mean.dir <- circ.mean(data)# Compute the mean direction of a random sample of observations. data <- runif(50, 0, pi) mean.dir <- circ.mean(data)
Creates a plot of circular data points on the current graphics device. Data points are either plotted as points on the unit circle, or the range of the circle is divided into a specified number of bins, and points are stacked in the bins corresponding to the number of observations in each bin.
circ.plot(x, main="", pch=16, stack=FALSE, bins=0, cex=1, dotsep=40, shrink=1)circ.plot(x, main="", pch=16, stack=FALSE, bins=0, cex=1, dotsep=40, shrink=1)
x |
vector of observations to be plotted, measured in radians. |
main |
title of plot. |
pch |
point character to use. See help on par. |
stack |
logical flag: if TRUE, points are stacked on the perimeter of the circle. Otherwise, all points are plotted on the perimeter of the circle. Default is FALSE. |
bins |
if stack = TRUE, bins is the number of arcs to partition the circle with. |
cex |
point character size. See help on par. |
dotsep |
constant used to specify the distance between stacked points, if stack = TRUE. Default is 40; larger values will create smaller spaces, while smaller values create larger spaces. This option can be useful when pch is other than 1, or when shrink is other than 1. |
shrink |
parameter that controls the size of the plotted circle. Default is 1. Larger values shrink the circle, while smaller values enlarge the circle. This option is useful when stack = TRUE. |
When there are many closely distributed observations, stacking is recommended. Otherwise, much information can be lost due to overplotting of points. When stacking the points, if there are many points in a particular bin, it may be necessary to shrink the plot of the circle so that all points fit. This is controlled with the parameter shrink. Generally the parameter dotsep does not need adjustment, however, when shrinking the plot, or for a very large number of observations, it may be helpful.
this function use function eqscplot from package MASS.
# Generate 100 observations from a von Mises distribution. # with mean direction 0 and concentration 3. data.vm <- rvm(100, 0, 3) # Plot data set. All points do not fit on plot. circ.plot(data.vm, stack=TRUE, bins=150) # Shrink the plot so that all points fit. circ.plot(data.vm, stack=TRUE, bins=150, shrink=1.5)# Generate 100 observations from a von Mises distribution. # with mean direction 0 and concentration 3. data.vm <- rvm(100, 0, 3) # Plot data set. All points do not fit on plot. circ.plot(data.vm, stack=TRUE, bins=150) # Shrink the plot so that all points fit. circ.plot(data.vm, stack=TRUE, bins=150, shrink=1.5)
Computes the circular range of a data set and performs a test of uniformity if specified.
circ.range(x, test=FALSE)circ.range(x, test=FALSE)
x |
vector of directional data measured in radians. |
test |
logical flag: if TRUE then the test of uniformity is performed; otherwise the test is not performed. Default is FALSE. |
The circular range is the shortest arc on the circle containing the entire set of data. The p-value is computed using the exact distribution of the circular range under the hypothesis of uniformity.
Returns a dataframe with the circular range, r. If the significance test is requested the p-value of the test is returned as p.value.
nCk(n,k) evaluate combinatoric, it should not be called by the users.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.4, World Scientific Press, Singapore.
kuiper, rao.spacing, r.test, v0.test, watson.
data <- rvm(50, 0, 2) circ.range(data, test=TRUE) data <- runif(50, 0, 2*pi) circ.range(data, test=TRUE)data <- rvm(50, 0, 2) circ.range(data, test=TRUE) data <- runif(50, 0, 2*pi) circ.range(data, test=TRUE)
Fits a regression model for a circular dependent and circular independent variable.
circ.reg(alpha, theta, order=1, level=0.05)circ.reg(alpha, theta, order=1, level=0.05)
alpha |
vector of data for the independent circular variable, measured in radians. |
theta |
vector of data for the dependent circular variable, measured in radians. |
order |
order of trigonometric polynomial to be fit. order must be an integer value. By default, order=1. |
level |
level of the test for the significance of higher order trigonometric terms. |
A trigonometric polynomial of alpha is fit against the cosine and sine of theta. The order of trigonometric polynomial is specified by order. Fitted values of theta are obtained by taking the inverse tangent of the predicted values of sin(theta) devided by the predicted values of cos(theta). Details of the regression model can be found in Sarma and Jammalamadaka (1993).
A data frame is returned with the following components.
rho |
square root of the average of the squares of the estimated conditional concentration parameters of theta given alpha. |
fitted |
fitted values of the model. |
data |
matrix whose columns correspond to alpha and theta. |
residuals |
circular residuals of the model. |
coef |
matrix whose entries are the estimated coefficients of the model. The first column corresponds to the coefficients of the model predicting the cosine of theta, while the second column contains the estimates for the model predicting the sine of theta. The rows of the matrix correspond to the coefficients according to increasing trigonometric order. |
pvalues |
p-values testing whether the (order + 1) trigonometric terms are significantly different from zero. |
A.k |
is mean of the cosines of the circular residuals. |
kappa |
assuming the circular residuals come from a von Mises distribution, kappa is the MLE of the concentration parameter. |
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 8.3, World Scientific Press, Singapore.
Sarma, Y. and Jammalamadaka, S. (1993). Circular Regression. Statistical Science and Data Analysis, 109-128. Proceeding of the Thrid Pacific Area Statistical Conference. VSP: Utrecht, Netherlands.
# Generate a data set of dependent circular variables. data1 <- runif(50, 0, 2*pi) data2 <- atan2(0.15*cos(data1) + 0.25*sin(data1), 0.35*sin(data1)) + rvm(50, 0, 5) # Fit a circular regression model. circ.lm <- circ.reg(data1, data2, order=1) # Obtain a crude plot a data and fitted regression line. plot(data1, data2) circ.lm$fitted[circ.lm$fitted>pi] <- circ.lm$fitted[circ.lm$fitted>pi] - 2*pi points(data1[order(data1)], circ.lm$fitted[order(data1)], type='l')# Generate a data set of dependent circular variables. data1 <- runif(50, 0, 2*pi) data2 <- atan2(0.15*cos(data1) + 0.25*sin(data1), 0.35*sin(data1)) + rvm(50, 0, 5) # Fit a circular regression model. circ.lm <- circ.reg(data1, data2, order=1) # Obtain a crude plot a data and fitted regression line. plot(data1, data2) circ.lm$fitted[circ.lm$fitted>pi] <- circ.lm$fitted[circ.lm$fitted>pi] - 2*pi points(data1[order(data1)], circ.lm$fitted[order(data1)], type='l')
Computes circular summary statistics including the sample size, mean direction and mean resultant length.
circ.summary(x)circ.summary(x)
x |
vector of data points measured in radians. |
Each observation is treated as a unit vector or a point on the unit circle. The resultant vector of the observations is found, and the direction of the resultant vector is returned as well as its length divided by the sample size.
Returns a data frame with variables n, the sample size; mean.dir, the sample mean direction; and rho, the sample mean resultant length.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 1.3, World Scientific Press, Singapore.
circ.mean, circ.disp, est.kappa, est.rho.
# Compute summary statistics of a random sample of observations. data <- runif(50, 0, pi) circ.summary(data)# Compute summary statistics of a random sample of observations. data <- runif(50, 0, pi) circ.summary(data)
Returns the cardoid density function evaluated at a particular value.
dcard(theta, mu, r)dcard(theta, mu, r)
theta |
vector of angles measured in radians at which to evaluate the density function. |
mu |
mean direction of the distribution. |
r |
concentration parameter of the distribution. Absolute value of r must be less than 0.5. |
Returns the cardoid density function evaluated at theta.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.2, World Scientific Press, Singapore.
Evaluates the density function of a mixture of two von Mises distributions.
dmixedvm(theta, mu1, mu2, kappa1, kappa2, p)dmixedvm(theta, mu1, mu2, kappa1, kappa2, p)
theta |
vector of values at which to evaluate the density function. |
mu1 |
mean direction in radians of one of the two von Mises distributions. |
mu2 |
mean direction in radians of the other von Mises distribution. |
kappa1 |
concentration parameter of one of the two von Mises distributions. |
kappa2 |
concentration parameter of the other von Mises distribution. |
p |
mixing proportion. |
Evaluates the density function , where VM is the von Mises density function.
Evaluates the density function of a mixture of two von Mises distributions at a given vector of values, theta.
Returns the triangular density function evaluated at a particular value.
dtri(theta, r)dtri(theta, r)
theta |
vector of angles measured in radians at which to evaluate the density function. |
r |
concentration parameter of the distribution. r must be between 0 and |
Returns the triangular density function evaluated at theta.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.3, World Scientific Press, Singapore.
Returns the von Mises density function evaluated at a particular value.
dvm(theta, mu, kappa)dvm(theta, mu, kappa)
theta |
vector of angles measured in radians at which to evaluate the density function. |
mu |
mean direction of the distribution measured in radians. |
kappa |
non-negative numeric value for the concentration parameter of the distribution. |
Returns the von Mises density function evaluated at theta.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.4, World Scientific Press, Singapore.
Returns the wrapped Cauchy density function evaluated at a particular value.
dwrpcauchy(theta, mu, rho)dwrpcauchy(theta, mu, rho)
theta |
vector of angles measured in radians at which to evaluate the density function. |
mu |
mean direction of the distribution measured in radians. |
rho |
concentration parameter of the distribution. rho must be in the interval from 0 to 1. |
Returns the wrapped Cauchy density function evaluated at theta.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.7, World Scientific Press, Singapore.
Estimate of the wrapped normal density function.
dwrpnorm(theta, mu, rho, sd=1, acc=1e-5, tol=acc)dwrpnorm(theta, mu, rho, sd=1, acc=1e-5, tol=acc)
theta |
value at which to evaluate the density function, measured in radians. |
mu |
mean direction of distribution, measured in radians. |
rho |
mean resultant length of distribution. |
sd |
different way of select |
acc |
parameter defining the accuracy of the estimation of the
density. Terms are added to the infinite summation that defines the
density function until successive estimates are within |
tol |
the same as |
The form of the wrapped normal density function is an infinite series
with index going from negative infinity to positive infinity. This
function begins with the zeroth term and adds terms to the series,
corresponding to both the positive and negative index, until the
summation changes by less than the parameter value of acc. You
can set rho by using sd with the following equivalence:
Returns an estimate of the wrapped normal density function.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.6, World Scientific Press, Singapore.
# Values for which to evaluate density theta <- c(1:500)*2*pi/500 #Compute wrapped normal density function density <- c(1:500) for(i in 1:500) density[i] <- dwrpnorm(theta[i], pi, .75) plot(theta, density) #Approximate area under density curve sum(density*2*pi/500)# Values for which to evaluate density theta <- c(1:500)*2*pi/500 #Compute wrapped normal density function density <- c(1:500) for(i in 1:500) density[i] <- dwrpnorm(theta[i], pi, .75) plot(theta, density) #Approximate area under density curve sum(density*2*pi/500)
Computes the maximum likelihood estimate of kappa, the concentration parameter of a von Mises distribution, given a set of angular measurements.
est.kappa(x, bias=FALSE)est.kappa(x, bias=FALSE)
x |
vector of circular data measured in radians. |
bias |
logical parameter determining whether a bias correction is used in the computation of the MLE. Default for bias is FALSE - the bias correction is not used. |
Best and Fisher (1981) show that the MLE of kappa is seriously biased when both sample size and mean resultant length are small. They suggest a bias-corrected estimate for kappa when n < 16.
Maximum likelihood estimate of kappa, the concentration parameter of a von Mises distribution, given a set of angular measurements.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 4.2.1, World Scientific Press, Singapore.
Best, D. and Fisher N. (1981). The bias of the maximum likelihood estimators of the von Mises-Fisher concentration parameters. Communications in Statistics - Simulation and Computation, B10(5), 493-502.
circ.mean, circ.disp, circ.summary, est.rho
data <- rvm(15, 0, 3) est.kappa(data) est.kappa(data, bias=TRUE)data <- rvm(15, 0, 3) est.kappa(data) est.kappa(data, bias=TRUE)
Returns the mean resultant length of a vector of circular data.
est.rho(x)est.rho(x)
x |
vector of data points measured in radians. |
Each observation is treated as a unit vector, or point on the unit circle. The resultant vector of the observations is found, and the length of the resultant vector divided by the sample size is returned.
Returns the mean resultant length of data.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 1.3, World Scientific Press, Singapore.
circ.mean, circ.disp, circ.summary, est.kappa.
# Compute the mean resultant length of a random sample of observations. data <- runif(100, 0, 2*pi) est.rho(data)# Compute the mean resultant length of a random sample of observations. data <- runif(100, 0, 2*pi) est.rho(data)
An alias of besselI(x, nu=0, expon.scaled = FALSE) used for compatibility with old version package code.
I.0(x)I.0(x)
x |
non-negative numerical value at which to evaluate the Bessel function. |
The use of this function is deprecated. Please use besselI.
Returns the zeroth order Bessel function of the first kind evaluated at a specified real number.
An alias of besselI(x, nu=1, expon.scaled = FALSE) used for compatibility with old version package code.
I.1(x)I.1(x)
x |
non-negative numerical value at which to evaluate the Bessel function. |
The use of this function is deprecated. Please use besselI.
Returns the first order Bessel function of the first kind, evaluated at a specified real number.
An alias of besselI(x, nu=p, expon.scaled = FALSE) used for compatibility with old version package code.
I.p(p, x)I.p(p, x)
p |
positive integer order of the Bessel function. |
x |
non-negative numerical value at which to evaluate the Bessel function. |
The use of this function is deprecated. Please use besselI.
Returns the p-th order Bessel function of the first kind, evaluated at a specified real number.
Performs Kuiper's one-sample test of uniformity on the circle.
kuiper(x, alpha=0)kuiper(x, alpha=0)
x |
vector of angular measurements in radians. |
alpha |
significance level of the test. Possible levels are 0.15, 0.1, 0.05, 0.025, 0.01. Alpha may be omitted or set to zero, in which case a range for the p-value of the test will be returned. |
Kuiper's test statistic is a rotation-invariant Kolmogorov-type test statistic. The critical values of a modified Kuiper's test statistic are used according to the tabulation given in Stephens (1970).
NULL
Kuiper's one-sample test of uniformity is performed, and the results are printed to the screen. If alpha is specified and non-zero, the test statistic is printed along with the critical value and decision. If alpha is omitted, the test statistic is printed and a range for the p-value of the test is given.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.2, World Scientific Press, Singapore.
Stephens, M. (1970). Use of the Kolmogorov-Smirnov, Cramer-von Mises and related statistics without extensive tables. Journal of the Royal Statistical Society, B32, 115-122.
circ.range, rao.spacing, r.test, v0.test, watson
# Generate data from the uniform distribution on the circle. data <- runif(100, 0, 2*pi) kuiper(data) # Generate data from the von Mises distribution. data <- rvm(100, 0, 3) kuiper(data, alpha=0.01)# Generate data from the uniform distribution on the circle. data <- runif(100, 0, 2*pi) kuiper(data) # Generate data from the von Mises distribution. data <- rvm(100, 0, 3) kuiper(data, alpha=0.01)
Plots the empirical distribution function of a data set.
plotedf(x, ...)plotedf(x, ...)
x |
vector of circular data measured in radians. |
... |
optional graphical parameters. See help section on par. |
The vector data is taken modulo 2*pi, and then the linear ranks are used to generate an empirical distribution function.
Creates a plot of the empirical distribution function of the vector data.
# Compare the edf's of two simulated sets of data. data1 <- rvm(10, 0, 3) data2 <- rvm(10, 0, 1) plotedf(data1, xlab="Data", ylab="EDF", main="Plots of Two EDF's") par(new=TRUE) plotedf(data2, axes=FALSE, xlab="", ylab="", lty=2)# Compare the edf's of two simulated sets of data. data1 <- rvm(10, 0, 3) data2 <- rvm(10, 0, 1) plotedf(data1, xlab="Data", ylab="EDF", main="Plots of Two EDF's") par(new=TRUE) plotedf(data2, axes=FALSE, xlab="", ylab="", lty=2)
Plots the empirical distribution of a data set against the best fitting von Mises distribution function.
pp.plot(x, ref.line=TRUE)pp.plot(x, ref.line=TRUE)
x |
vector of angular measurements in radians. |
ref.line |
logical flag: if TRUE a 45 degree reference line is added to the plot. Default is TRUE. |
The maximum likelihood estimates of the parameters of the von Mises distribution are computed from the given data set. The empirical distribution function is plotted against a von Mises distribution function with parameters given by the MLEs computed.
NULL
A plot is created on the current graphics device.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 10.2, World Scientific Press, Singapore.
data <- rvm(25, 0, 3) pp.plot(data) data <- c(rvm(20, 0, 7), rvm(20, pi, 7)) pp.plot(data)data <- rvm(25, 0, 3) pp.plot(data) data <- c(rvm(20, 0, 7), rvm(20, pi, 7)) pp.plot(data)
Estimates the cummulative probability for a von Mises distribution.
pvm(theta, mu, kappa, acc=1e-020)pvm(theta, mu, kappa, acc=1e-020)
theta |
angular value in radians. |
mu |
mean direction of the von Mises distribution. |
kappa |
concentration parameter of the von Mises distribution. |
acc |
parameter relating to the accuracy of the estimated cummulative probability. See details below. Default value is 1e-020. |
Cummulative probabilities are computed according to the expression for the von Mises cdf given by Gumbel et al. (1953), which gives the cdf as a function of an infinite sum. The parameter acc specifies the accuracy with which this sum is approximated. Terms greater than acc are included in the summation. The cummulative probabilities given by pvm coincide with those tabulated by Mardia (1972), which are given to five decimal values.
Returns the probability that a von Mises random variable falls between 0 and theta.
Gumbel, E., Greenwood, J. and Durand, D. (1953). The circular normal distribution: theory and tables. Journal of the American Statistical Association, 48, 131-152.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.4, World Scientific Press, Singapore.
Performs a Rayleigh test of uniformity, assessing the significance of the mean resultant length. The alternative hypothesis is a unimodal distribution with unknown mean direction and unknown mean resultant length.
r.test(x, degree=FALSE)r.test(x, degree=FALSE)
x |
vector of angular measurements in radians. |
degree |
logical flag: if TRUE, data is assumed to be measured in degrees rather than radians. Default is FALSE. |
Returns a list with two components: the mean resultant length, r.bar, and the p-value of the test statistic, p.value.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Sections 3.3.2 and 3.4.1, World Scientific Press, Singapore.
circ.range, kuiper, rao.spacing, v0.test, watson
data <- rvm(25, pi, 2) r.test(data)data <- rvm(25, pi, 2) r.test(data)
Performs Rao's test for homogeneity on k populations of angular data.
rao.homogeneity(x, alpha=0)rao.homogeneity(x, alpha=0)
x |
a list of angular variables measured in radians. Each component of the list should be a vector corresponding to one of the k samples. The data must be in the form of a list, to allow for samples of varying sizes. |
alpha |
numeric value specifying the significance level of the test. Default is 0, in which case p-values for the test statistic are returned. |
Critical values and p-values are determined according to the chi-squared approximation of the the test statistic.
NULL
The test is performed, and the results are written to the screen. Test results are given for both the test of equality of polar vectors, and of dispersions. If alpha is specified, the test statistic is returned, along with the level alpha critical value. If alpha is not specified, a p-value for the test is computed.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.6.1, World Scientific Press, Singapore.
Rao, J.S. (1967). Large sample tests for the homogeneity of angular data, Sankhya, Ser, B., 28.
Performs Rao's Spacing Test of uniformity.
rao.spacing(x, alpha=0, rad=TRUE)rao.spacing(x, alpha=0, rad=TRUE)
x |
numeric vector of angular data measured in degrees. |
alpha |
numeric value specifying the significance level of the test. The default value is 0, in which case, a range for the p-value will be returned. Valid significance levels are 0.10, 0.05, 0.01 and 0.001. |
rad |
logical value. If true, data will be assumed to be measured radians. If false, data will be assumed to be measured in degrees. |
If alpha is specified, critical values are determined from a table of simulated critical points (see reference below). If alpha is not specified, a range for the p-value is determined using the table of simulated critical points.
NULL
The Rao's Spacing Test is performed, and the results are written to the screen.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.4, World Scientific Press, Singapore.
Rao, J.S. (1976). Some tests based on arc-lengths for the circle. Sankhya, The Indian Journal of Statistics, Serial B(4), 38, 329-338.
Russell, G.S. and Levitin, D.J. (1995). An expanded table of probability values for Rao's Spacing Test. Communications in Statistics - Simulation and Computation, 24, 4, 879-888.
circ.range, kuiper, r.test, v0.test, watson
Table for Rao spacing function
data(rao.table)data(rao.table)
Generates pseudo-random numbers from a cardoid distribution.
rcard(n, mu, r)rcard(n, mu, r)
n |
number of random variables to generate. |
mu |
mean direction of the distribution in radians. |
r |
concentration parameter of the distribution. The absolute value of r must be less than 0.5. |
An acceptance/rejection method of simulation is employed.
Returns a vector of n independent random variables generated from a cardoid distribution.
Generates pseudo-random numbers from a mixture of two von Mises distributions.
rmixedvm(n, mu1, mu2, kappa1, kappa2, p)rmixedvm(n, mu1, mu2, kappa1, kappa2, p)
n |
number of random variables to generate. |
mu1 |
mean direction in radians of one of the two von Mises distributions. |
mu2 |
mean direction in radians of the other von Mises distribution. |
kappa1 |
concentration parameter of one of the two von Mises distributions. |
kappa2 |
concentration parameter of the other von Mises distribution. |
p |
mixing proportion. |
Simulates random variables with the density function p*VM(mu1, kappa1) + (1-p)*VM(mu2, kappa2), where VM is the von Mises density function.
A uniform random variable on (0,1) is generated. If it is less than p, then a variable is generated from VM(mu1, kappa1), else a variable is generated from VM(mu2, kappa2). Simulation from the von Mises distribution is done via the algorithm due to Best and Fisher (1979).
Returns a vector of n independent random variables generated from a mixture of two von Mises distributions.
Best, D. and Fisher, N. (1979). Efficient simulation of the von Mises distribution. Applied Statistics, 24, 152-157.
Creates a rose diagram of a circular data set on the current graphics device.
rose.diag(x, bins, main="", prop=1, pts=FALSE, cex=1, pch=16, dotsep=40, shrink=1)rose.diag(x, bins, main="", prop=1, pts=FALSE, cex=1, pch=16, dotsep=40, shrink=1)
x |
vector of directional data measured in radians. |
bins |
number of groups to partition the circle with. This will be the number of petals or sectors in the rose diagram. |
main |
title of plot. Default is no title. |
prop |
numerical constant determining the radii of the sectors. By default, prop = 1. |
pts |
logical value. If TRUE, points will be stacked on the circumference of the circle, according to bins - one point for each observation. The default value is FALSE, no points plotted. |
cex |
size of points, if pts = TRUE. By default, cex = 1. See help on cex in help section for par. |
pch |
if pts = TRUE, pch determines the character used. See help on pch in help section for par. |
dotsep |
constant used to specify the distance between stacked points, if pts = TRUE. Default is 40; larger values will create smaller spaces, while smaller values create larger spaces. This option can be useful when pch is other than 1, or when shrink is other than 1 (see below). |
shrink |
parameter that controls the size of the plotted circle. Default is 1. Larger values shrink the circle, while smaller values enlarge the circle. This option can be useful when pts = TRUE. |
The circumference of the circle is split into groups, the number of groups specified by bins. For each group, a sector is drawn. The radii of the sectors are by default equal to the square root of the relative frequencies of observations in each group. This ensures that the area of the sector is proportional to the group frequency. The length of the radii can be controlled by varying the parameter prop.
NULL
A rose diagram is plotted on the current graphics device.
# Generate uniform data and create several rose diagrams. # Some optional parameters may be needed to optimize plots. data <- runif(50, 0, 2*pi) rose.diag(data, bins = 18, main = 'Uniform Data') rose.diag(data, bins = 18, main = 'Stacked Points', pts=TRUE) # Generate von Mises data and create several rose diagrams. data <- rvm(25, 0, 5) rose.diag(data, bins=18, pts=TRUE) # Points fall out of bounds. rose.diag(data, bins=36, prop=1.5, pts=TRUE, shrink=1.5) # Adjust optional parameters to fit all points on plot.# Generate uniform data and create several rose diagrams. # Some optional parameters may be needed to optimize plots. data <- runif(50, 0, 2*pi) rose.diag(data, bins = 18, main = 'Uniform Data') rose.diag(data, bins = 18, main = 'Stacked Points', pts=TRUE) # Generate von Mises data and create several rose diagrams. data <- rvm(25, 0, 5) rose.diag(data, bins=18, pts=TRUE) # Points fall out of bounds. rose.diag(data, bins=36, prop=1.5, pts=TRUE, shrink=1.5) # Adjust optional parameters to fit all points on plot.
Returns random deviates from the stable family of probability distributions.
rstable(n, scale = 1, index = stop("no index arg"), skewness = 0)rstable(n, scale = 1, index = stop("no index arg"), skewness = 0)
n |
sample size |
index |
number from the interval (0, 2]. An index of 2 corresponds to the normal, 1 to the Cauchy. Smaller values mean longer tails. |
skewness |
number giving the modified skewness (see Chambers et al., 1976). Negative values correspond to skewness to the left (the median is smaller than the mean, if it exists), and positive values correspond to skewness to the right (the median is larger than the mean). The absolute value of skewness should not exceed 1. |
scale |
the scale of the distribution. |
This function returns a random variate from the Levy skew stable
distribution with index=, scale= and skewness=. The skewness parameter must lie in the range [-1,1] while the index parameter must lie in the range (0,2]. The Levy skew stable probability distribution is defined by a fourier transform,
When the term is replaced by . For the distribution reduces to a Gaussian distribution with and the skewness parameter has no effect. For the tails of the distribution become extremely
wide. The symmetric distribution corresponds to .
The Levy alpha-stable distributions have the property that if alpha-stable variates are drawn from the distribution then the sum $Y = X_1 + X_2 + ... + X_N$ will also be distributed as an alpha-stable variate, .
There is no explicit solution for the form of and there are no density, probability or quantile functions supplied for this distribution.
random sample from the specified stable distribution.
Chambers, J. M., Mallows, C. L. and Stuck, B. W. (1976). A Method for Simulating Stable Random Variables. Journal of the American Statistical Association 71, 340-344.
Lo\"gae\"ve, M. (1977). Probability Theory I. (fourth edition) Springer-Verlag, New York.
hist(rstable(200, 1.5, .5)) #fairly long tails, skewed righthist(rstable(200, 1.5, .5)) #fairly long tails, skewed right
Generates pseudo-random numbers from a symmetric triangular distribution.
rtri(n, r)rtri(n, r)
n |
number of random variables to generate. |
r |
concentration parameter of the distribution. r must be between 0 and |
The triangular cdf is inverted to obtain the random numbers.
Returns a vector of n independent random variables generated from a symmetric triangular distribution with mean direction 0 and concentration parameter r.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.3.3, World Scientific Press, Singapore.
Generates pseudo-random numbers from a von Mises distribution.
rvm(n, mean, k)rvm(n, mean, k)
n |
number of random variables to generate. |
mean |
mean direction in radians of the von Mises distribution. |
k |
concentration parameter of the von Mises distribution. |
Simulation from the von Mises distribution is done via the algorithm due to Best and Fisher (1979).
Returns a vector of n independent random variables generated from a von Mises distribution. Values are between 0 and 2 pi.
Best, D. and Fisher, N. (1979). Efficient simulation of the von Mises distribution. Applied Statistics, 24, 152-157.
Generates pseudo-random numbers from a wrapped Cauchy distribution.
rwrpcauchy(n, location=0, rho=exp(-1))rwrpcauchy(n, location=0, rho=exp(-1))
n |
number of random variables to generate. |
location |
modal angle in radians of the wrapped Cauchy distribution. |
rho |
concentration parameter. rho must be between 0 and 1. |
n random variables are generated from a Cauchy distribution with location parameter location, and concentration parameter -log(rho). The function returns these values modulo 2*pi.
Returns a vector of n independent random variables generated from a wrapped Cauchy distribution.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.6, World Scientific Press, Singapore.
Generates pseudo-random numbers from a wrapped normal distribution.
rwrpnorm(n, mu, rho, sd=1)rwrpnorm(n, mu, rho, sd=1)
n |
number of random variables to generate. |
mu |
mean direction in radians of the wrapped normal distribution. |
rho |
concentration parameter. rho must be between 0 and 1. |
sd |
different way of select |
n random variables are generated from a normal distribution with mean
direction mu, and variance -2*log(rho). The function returns these
values modulo 2*pi. You can set rho by using sd with the following equivalence:
Returns a vector of n independent random variables generated from a wrapped normal distribution.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.5, World Scientific Press, Singapore.
Generates pseudo-random numbers from a wrapped stable distribution.
rwrpstab(n, index, skewness, scale=1)rwrpstab(n, index, skewness, scale=1)
n |
number of random variables to generate. |
index |
number from the interval (0, 2]. An index of 2 corresponds to the normal, 1 to the Cauchy. Smaller values mean longer tails. |
skewness |
number giving the modified skewness. Negative values correspond to skewness to the left (the median is smaller than the mean, if it exists), and positive values correspond to skewness to the right (the median is larger than the mean). The absolute value of skewness should not exceed 1. |
scale |
the scale of the distribution |
n random variables are generated from a stable distribution with with parameters index, skewness and scale. The function returns these values modulo 2*pi.
Returns a vector of n independent random variables generated from a wrapped stable distribution.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 2.2.8, World Scientific Press, Singapore.
Computes the specified order trigonometric moment for a set of directional data points.
trig.moment(x, p=1, center=FALSE)trig.moment(x, p=1, center=FALSE)
x |
directional data set measured in radians. |
p |
order of trigonometric moment to be computed. Default is for a first order trigonometric moment. |
center |
logical flag whether to compute centered moments or not. Default is to compute an uncentered moment. |
Returns a data frame with variables mu.p, rho.p, cos.p, and sin.p, respectively the pth trigonometric moment's direction, resultant length, and real and imaginary components.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 1.3, World Scientific Press, Singapore.
circ.disp, circ.mean, circ.summary, est.kappa, est.rho
Performs a Rayleigh test of uniformity, assessing the significance of the mean resultant length. The alternative hypothesis is a unimodal distribution with a specified mean direction and unknown mean resultant length.
v0.test(x, mu0=0, degree=FALSE)v0.test(x, mu0=0, degree=FALSE)
x |
vector of angular measurements in radians. |
mu0 |
Specified mean direction in alternative hypothesis, measured in the same units (degrees or radians) as the data. |
degree |
logical flag: if TRUE, data is assumed to be measured in degrees rather than radians. Default is FALSE. |
Returns a list with two components: the mean resultant length, r0.bar, and the p-value of the test statistic, p.value.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Sections 3.3.3 and 3.4.1, World Scientific Press, Singapore.
circ.range, kuiper, rao.spacing, r.test, watson
data <- rvm(25, pi, 2) v0.test(data, mu0=pi)data <- rvm(25, pi, 2) v0.test(data, mu0=pi)
Generates simple bootstrap confidence intervals for the parameters of a von Mises distribtution: the mean direction mu, and the concentration parameter kappa.
vm.bootstrap.ci(x, bias=FALSE, alpha=0.05, reps=1000, print=TRUE)vm.bootstrap.ci(x, bias=FALSE, alpha=0.05, reps=1000, print=TRUE)
x |
vector of angular measurements in radians. |
bias |
logical flag: if TRUE, the replication estimates for kappa are computed with a bias corrected method. See est.kappa. Default is FALSE, i.e. no bias correction. |
alpha |
parameter determining level of confidence intervals. 1-alpha confidence intervals for mu and kappa are computed. By default, 95% confidence intervals are generated. |
reps |
number of resampled data sets to use. Default is 1000. |
print |
logical flag indicating whether the algorithm should print a message indicating which set of replicates is currently being drawn. Default is TRUE. |
Percentile confidence intervals are computed by resampling from the original data set B times. For each resampled data set, the MLE's of mu and kappa are computed. The bootstrap confidence intervals are the alpha/2 and 1-alpha/2 percentiles of the B MLE's computed for each resampled data set.
A list is returned with the following components: mu.ci and kappa.ci contain the limits of the confidence intervals for mu and kappa respectively. mu.reps and kappa.reps contain the estimates of mu and kappa for each resampled data set.
The confidence intervals are printed to the screen.
x <- rvm(25, 0, 3) x.bs <- vm.bootstrap.ci(x, alpha=.10) hist(x.bs$kappa.reps)x <- rvm(25, 0, 3) x.bs <- vm.bootstrap.ci(x, alpha=.10) hist(x.bs$kappa.reps)
Computes the maximum likelihood estimates for the parameters of a von Mises distribution: the mean direction and the concentration parameter.
vm.ml(x, bias=FALSE)vm.ml(x, bias=FALSE)
x |
vector of angular measurements in radians. |
bias |
logical flag: if TRUE, the estimate for kappa is computed with a bias corrected method. See est.kappa. Default is FALSE, i.e. no bias correction. |
Returns a data frame with two components: mu and kappa, the MLES of the respective parameters.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 4.2.1, World Scientific Press, Singapore.
circ.mean, est.kappa, vm.bootstrap.ci
Performs a Watson's goodness of fit test for the von Mises or circular uniform distribution.
watson(x, alpha=0, dist='uniform')watson(x, alpha=0, dist='uniform')
x |
vector of angular measurements in radians. |
alpha |
significance level of the test. Valid levels are 0.01, 0.05, 0.1. This argument may be ommited, in which case, a range for the p-value will be returned. |
dist |
distribution to test for. The default is the uniform distribution. To test for the von Mises distribution, set dist = 'vm'. |
NULL
If dist = 'uniform', Watson's one-sample test for the circular uniform distribution is performed, and the results are printed to the screen. If alpha is specified and non-zero, the test statistic is printed along with the critical value and decision. If alpha is omitted, the test statistic is printed and a range for the p-value of the test is given.
If dist = 'vm', estimates of the population parameters are used to evaluate the von Mises distribution function at all data points, thereby arriving at a sample of approximately uniformly distributed data, if the original observations have a von Mises distribution. The one-sample Watson test is then applied to the transformed data as above.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.2, World Scientific Press, Singapore.
Stephens, M. (1970). Use of the Kolmogorov-Smirnov, Cramer-von Mises and related statistics without extensive tables. Journal of the Royal Statistical Society, B32, 115-122.
circ.range, kuiper, rao.spacing, r.test, v0.test
# Generate data from the uniform distribution on the circle. data <- runif(100, 0, 2*pi) watson(data) # Generate data from a von Mises distribution. data <- rvm(50, 0, 4) watson(data, 0.05, dist='vm')# Generate data from the uniform distribution on the circle. data <- runif(100, 0, 2*pi) watson(data) # Generate data from a von Mises distribution. data <- rvm(50, 0, 4) watson(data, 0.05, dist='vm')
Performs Watson's test for homogeneity on two samples of circular data.
watson.two(x, y, alpha=0, plot=FALSE)watson.two(x, y, alpha=0, plot=FALSE)
x |
vector of circular data measured in radians. |
y |
vector of circular data measured in radians. |
alpha |
significance level of the test. Valid levels are 0.001, 0.01, 0.05, 0.1. This argument may be ommited, in which case, a range for the p-value will be returned. |
plot |
logical value. If TRUE, an overlayed plot of both empirical distribution functions will be sent to the current graphics device. The default value if FALSE. |
Critical values for the test statistic are obtained using the asymptotic distribution of the test statistic. It is recommended to use the obtained critical values and ranges for p-values only for combined sample sizes in excess of 17. Tables are available for smaller sample sizes and can be found in Mardia (1972) for instance.
NULL
Watson's two-sample test of homogeneity is performed, and the results are printed to the screen. If alpha is specified and non-zero, the test statistic is printed along with the critical value and decision. If alpha is omitted, the test statistic is printed and a range for the p-value of the test is given. If plot=TRUE, an overlayed plot of both empirical distribution functions will be sent to the current graphics device.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 7.5, World Scientific Press, Singapore.
# Perform a two-sample test of homogeneity on two # simulated data sets. data1 <- rvm(20, 0, 3) data2 <- rvm(20, pi, 2) watson.two(data1, data2, alpha=0.05, plot=TRUE) watson.two(data1, data2)# Perform a two-sample test of homogeneity on two # simulated data sets. data1 <- rvm(20, 0, 3) data2 <- rvm(20, pi, 2) watson.two(data1, data2, alpha=0.05, plot=TRUE) watson.two(data1, data2)
Computes the maximum likelihood estimates of the location and scale parameters for a wrapped Cauchy distribution.
wrpcauchy.ml(x, mu0, rho0, acc=1e-015)wrpcauchy.ml(x, mu0, rho0, acc=1e-015)
x |
vector of angular data measured in radians. |
mu0 |
initial estimate of the location parameter. |
rho0 |
initial estimate of the scale parameter. rho0 should be between in [0,1). |
acc |
degree of accuracy in the approximation of the estimates. The default value is 1e-15. See below for details. |
An iterative algorithm due to Kent and Tyler (1988) is used. Initial values of the MLE's are required. In the estimation process, estimates of quantities mu1 and mu2 (see Kent and Tyler) are updated iteratively. When both mu1 and mu2 change by less than acc from one iteration to the next, the process ends, and mu and rho are computed from the final estimates of mu1 and mu2.
A dataframe is returned with the variables mu and rho, where mu and rho are the approximations of the MLE of the location and scale parameter, respectively.
Kent, J. and Tyler, D. (1988). Maximum likelihood estimation for the wrapped Cauchy distribution. Journal of Applied Statistics, 15, 2, 247-254.
Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 4.4, World Scientific Press, Singapore.
# Generate data from a wrapped Cauchy distribution. data <- rwrpcauchy(50, 0, 0.75) # Compute the sample mean direction and resultant length. mu0 <- circ.mean(data) rho0 <- est.rho(data) # Estimate the MLE's of the Cauchy distribution parameters. wrpcauchy.ml(data, mu0, rho0)# Generate data from a wrapped Cauchy distribution. data <- rwrpcauchy(50, 0, 0.75) # Compute the sample mean direction and resultant length. mu0 <- circ.mean(data) rho0 <- est.rho(data) # Estimate the MLE's of the Cauchy distribution parameters. wrpcauchy.ml(data, mu0, rho0)