Geostatistics with R: Directional Statistical Preliminaries with R
Dr. Debashis Chatterjee
5/29/2024
# Install the required packages
#install.packages("circular")
#install.packages("ggplot2")# Load the libraries
library(circular)##
## Attaching package: 'circular'## The following objects are masked from 'package:stats':
##
## sd, varlibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3Simulate a Geological Dataset
Let’s create a synthetic dataset with directional data such as strike and dip measurements.
# Set seed for reproducibility
set.seed(123)
# Simulate directional data
n <- 100
strike <- circular(runif(n, 0, 360), units = "degrees")
dip <- runif(n, 0, 90)
# Convert dip data to a circular object
dip_circular <- circular(dip, units = "degrees")
# Create a data frame
geo_data <- data.frame(
ID = 1:n,
Strike = strike,
Dip = dip,
dip_circular=dip_circular
)
# View the first few rows of the dataset
head(geo_data)## ID Strike Dip dip_circular
## 1 1 103.52791 53.99901 53.99901
## 2 2 283.78985 29.95412 29.95412
## 3 3 147.23169 43.97517 43.97517
## 4 4 317.88627 85.90264 85.90264
## 5 5 338.56822 43.46122 43.46122
## 6 6 16.40034 80.13152 80.13152Visualize Strike and Dip Together
Create a scatter plot to visualize strike and dip together.
# Convert circular data to numeric for plotting
geo_data$Strike_num <- as.numeric(geo_data$Strike)
# Scatter plot of Strike vs Dip
ggplot(geo_data, aes(x = Strike_num, y = Dip)) +
geom_point(color = "blue", alpha = 0.6) +
labs(title = "Scatter Plot of Strike vs Dip", x = "Strike (degrees)", y = "Dip (degrees)") +
theme_minimal()
# Install the required packages
#install.packages("circular")
#install.packages("rgl")
# Load the libraries
library(circular)
library(rgl)## Warning: package 'rgl' was built under R version 4.3.3Convert Strike and Dip to Cartesian Coordinates
To plot on a sphere, convert the strike and dip to Cartesian coordinates. Assume the dip angle is from the horizontal plane.
# Function to convert strike and dip to Cartesian coordinates
convert_to_cartesian <- function(strike, dip) {
strike_rad <- strike * pi / 180
dip_rad <- dip * pi / 180
x <- cos(dip_rad) * cos(strike_rad)
y <- cos(dip_rad) * sin(strike_rad)
z <- sin(dip_rad)
return(data.frame(x, y, z))
}
# Convert the data
cartesian_coords <- convert_to_cartesian(geo_data$Strike, geo_data$Dip)
geo_data <- cbind(geo_data, cartesian_coords)
# View the first few rows of the dataset with Cartesian coordinates
head(geo_data)## ID Strike Dip dip_circular Strike_num x y
## 1 1 103.52791 53.99901 53.99901 103.52791 -0.13749739 0.57149144
## 2 2 283.78985 29.95412 29.95412 283.78985 0.20652240 -0.84145212
## 3 3 147.23169 43.97517 43.97517 147.23169 -0.60512151 0.38950065
## 4 4 317.88627 85.90264 85.90264 317.88627 0.05300373 -0.04791563
## 5 5 338.56822 43.46122 43.46122 338.56822 0.67565071 -0.26521700
## 6 6 16.40034 80.13152 80.13152 16.40034 0.16441379 0.04839067
## z
## 1 0.8090068
## 2 0.4993063
## 3 0.6943466
## 4 0.9974441
## 5 0.6878634
## 6 0.9852038# Open a new 3D plotting device
open3d()## wgl
## 1# Plot the sphere
spheres3d(0, 0, 0, radius = 1, color = "lightblue", alpha = 0.3, shininess = 50)
# Plot the data points
points3d(geo_data$x, geo_data$y, geo_data$z, color = "blue", size = 5)
# Add axes for better understanding
axes3d(c("x-", "y-", "z-"), col = "black")
title3d(xlab = "X", ylab = "Y", zlab = "Z", main = "3D Scatter Plot of Strike and Dip on Sphere")
# Add grid lines
grid3d(c("x", "y", "z"))Stereographic Projection
Convert Strike and Dip to Stereographic Projection Convert the strike and dip data to x and y coordinates using the stereographic projection.
# Function to convert strike and dip to stereographic projection
stereographic_projection <- function(strike, dip) {
rad_strike <- strike * pi / 180
rad_dip <- dip * pi / 180
x <- tan((pi / 2 - rad_dip) / 2) * sin(rad_strike)
y <- tan((pi / 2 - rad_dip) / 2) * cos(rad_strike)
return(data.frame(x = x, y = y))
}
# Apply the projection to the data
projection <- stereographic_projection(geo_data$Strike, geo_data$Dip)
geo_data$x <- projection$x
geo_data$y <- projection$y# Scatter plot of the stereographic projection
ggplot(geo_data, aes(x = x, y = y)) +
geom_point(color = "blue", alpha = 0.6) +
coord_fixed() +
labs(title = "Stereographic Projection of Strike and Dip", x = "X", y = "Y") +
theme_minimal()
Stereographic projection is a method of mapping points from a sphere onto a plane. This type of projection is conformal, meaning it preserves angles, which makes it particularly useful in fields such as geology for visualizing orientations and structural data.
Consider a sphere centered at the origin with radius \(R\). The north pole of the sphere is defined as the point \((0, 0, R)\). The stereographic projection maps points on the sphere (excluding the north pole) to points on the equatorial plane \(z = 0\).
Let \(P = (x, y, z)\) be a point on the sphere. The stereographic projection projects \(P\) onto the plane at a point \(P' = (X, Y, 0)\).
\end{figure}
To derive the coordinates \((X, Y)\) on the plane in terms of \((x, y, z)\) on the sphere, we start with the equation of the sphere: \[ x^2 + y^2 + z^2 = R^2 \]
Since \(P\) and \(P'\) lie on a line through the north pole, their coordinates are related by a scaling factor \(k\). Thus, we have: \[ P' = k P = k (x, y, z) \]
For \(P'\) to lie on the plane \(z = 0\), the z-coordinate of \(P'\) must be zero. Hence: \[ k z = 0 \]
Since \(P \neq\) North Pole, \(z \neq 0\), which implies \(k = 0\). However, the actual scaling factor involves using the line passing through the north pole: \[ P' = \left( \frac{x}{1 - \frac{z}{R}}, \frac{y}{1 - \frac{z}{R}}, 0 \right) \]
This transformation projects the point \(P\) on the sphere to the point \(P'\) on the plane.
The stereographic projection has several important properties:In geology, the stereographic projection is commonly used to plot the orientation of geological structures such as faults, bedding planes, and fold axes. The strike and dip of a plane can be represented as points on a stereographic projection, facilitating the analysis of spatial relationships and orientation patterns.
Stereographic projection is a powerful tool in geology and other fields that require the visualization of spherical data. Its conformal nature and ability to preserve circular structures make it particularly useful for representing orientations and structural features.
Visualize the Strike Data (Rose Diagram)
A rose diagram (circular histogram) is commonly used to visualize strike data.
# Rose diagram for strike data
rose_diag <- rose.diag(strike, bins = 36, prop = 1.5, col = "skyblue", main = "Rose Diagram of Strike Data")
Visualize the Dip Data
Dip data can be visualized using a histogram, but since it’s directional data, we use a circular plot.
# Circular plot for dip data
# Histogram for dip data
ggplot(geo_data, aes(x = Dip)) +
geom_histogram(binwidth = 5, fill = "coral", color = "black", alpha = 0.7) +
labs(title = "Histogram of Dip Data", x = "Dip (degrees)", y = "Frequency") +
theme_minimal()
# Create a rose diagram for dip data
rose_diag_dip <- rose.diag(geo_data$dip_circular, bins = 18, prop = 2, col = "coral")
Statistical Analysis of Strike Data
Perform a Rayleigh test to check if the strike data is uniformly distributed.
# Rayleigh test for uniformity
rayleigh_test <- rayleigh.test(strike)
# Print the test result
cat("Rayleigh Test Statistic:", rayleigh_test$statistic, "\n")## Rayleigh Test Statistic: 0.02338578cat("P-value:", rayleigh_test$p.value, "\n")## P-value: 0.9467791Circular Analog of the Normal Distribution:
The von Mises distribution is the most common probability distribution for circular data. It’s characterized by two parameters:
Mean Direction (\(\mu\)): The average or central direction of the data (in radians). Concentration Parameter (\(\kappa)\): A measure of how concentrated the data is around the mean direction. Higher values of kappa indicate more tightly clustered data.
Density Function:
# Simulate data from von Mises distribution
set.seed(123)
data <- rvonmises(n = 100, mu = circular(pi/4), kappa = 2) # Mean direction: 45 degrees, kappa = 2
# Plot the data
plot(data)
Circular-linear regression is used when one variable is circular (like strike or angle) and the other is linear (like dip or depth). In R, we can use the circular package for circular data and the circular package’s lm.circular function for performing circular-linear regression.
Here is the R code to perform circular-linear regression on the simulated dataset where strike is the circular variable and dip is the linear variable:
# Install the required packages
#install.packages("circular")
#install.packages("ggplot2")
#install.packages("circularRegression")
# Load the libraries
library(circular)
library(ggplot2)Circular-Linear Regression
Use the lm.circular function to perform the circular-linear regression.
# Generate a data set of dependent circular variables.
x <- circular(runif(50, 0, 2*pi))
y <- atan2(0.15*cos(x) + 0.25*sin(x), 0.35*sin(x)) +
rvonmises(n=50, mu=circular(0), kappa=5)
# Fit a circular-circular regression model.
circ.lm <- lm.circular(y, x, order=1)
# Obtain a crude plot of the data and fitted regression line.plot.default(x, y)
circ.lm$fitted[circ.lm$fitted>pi] <- circ.lm$fitted[circ.lm$fitted>pi] - 2*pi
points.default(x[order(x)], circ.lm$fitted[order(x)], type='l')
# Fit a circular-linear regression model and show predictions.
set.seed(1234)
x <- cbind(rnorm(10), rep(1,10))
y <- circular(2*atan(c(x%*%c(5,1))))+rvonmises(10, mu=circular(0), kappa=100)
lm.circular(y=y, x=x, init=c(5,1), type='c-l', verbose=TRUE)## Iteration 1 : Log-Likelihood = 24.56042
## Iteration 2 : Log-Likelihood = 24.56847
## Iteration 3 : Log-Likelihood = 24.57257
## Iteration 4 : Log-Likelihood = 24.57465
## Iteration 5 : Log-Likelihood = 24.5757
## Iteration 6 : Log-Likelihood = 24.57622
## Iteration 7 : Log-Likelihood = 24.57649
## Iteration 8 : Log-Likelihood = 24.57662
## Iteration 9 : Log-Likelihood = 24.57669
## Iteration 10 : Log-Likelihood = 24.57672
## Iteration 11 : Log-Likelihood = 24.57674
## Iteration 12 : Log-Likelihood = 24.57675
## Iteration 13 : Log-Likelihood = 24.57675
## Iteration 14 : Log-Likelihood = 24.57675
## Iteration 15 : Log-Likelihood = 24.57675
## Iteration 16 : Log-Likelihood = 24.57675
## Iteration 17 : Log-Likelihood = 24.57676
## Iteration 18 : Log-Likelihood = 24.57676
## Iteration 19 : Log-Likelihood = 24.57676
## Iteration 20 : Log-Likelihood = 24.57676
## Iteration 21 : Log-Likelihood = 24.57676
## Iteration 22 : Log-Likelihood = 24.57676
## Iteration 23 : Log-Likelihood = 24.57676
## Iteration 24 : Log-Likelihood = 24.57676
## Iteration 25 : Log-Likelihood = 24.57676
## Iteration 26 : Log-Likelihood = 24.57676
## Iteration 27 : Log-Likelihood = 24.57676
## Iteration 28 : Log-Likelihood = 24.57676
## Iteration 29 : Log-Likelihood = 24.57676
## Iteration 30 : Log-Likelihood = 24.57676
## Iteration 31 : Log-Likelihood = 24.57676
## Iteration 32 : Log-Likelihood = 24.57676
## Iteration 33 : Log-Likelihood = 24.57676
## Iteration 34 : Log-Likelihood = 24.57676
## Iteration 35 : Log-Likelihood = 24.57676
## Iteration 36 : Log-Likelihood = 24.57676
## Iteration 37 : Log-Likelihood = 24.57676
## Iteration 38 : Log-Likelihood = 24.57676
## Iteration 39 : Log-Likelihood = 24.57676
## Iteration 40 : Log-Likelihood = 24.57676
## Iteration 41 : Log-Likelihood = 24.57676
## Iteration 42 : Log-Likelihood = 24.57676
## Iteration 43 : Log-Likelihood = 24.57676
## Iteration 44 : Log-Likelihood = 24.57676
## Iteration 45 : Log-Likelihood = 24.57676
## Iteration 46 : Log-Likelihood = 24.57676
## Iteration 47 : Log-Likelihood = 24.57676
## Iteration 48 : Log-Likelihood = 24.57676
## Iteration 49 : Log-Likelihood = 24.57676
## Iteration 50 : Log-Likelihood = 24.57676
## Iteration 51 : Log-Likelihood = 24.57676
## Iteration 52 : Log-Likelihood = 24.57676
## Iteration 53 : Log-Likelihood = 24.57676
## Iteration 54 : Log-Likelihood = 24.57676
## Iteration 55 : Log-Likelihood = 24.57676
## Iteration 56 : Log-Likelihood = 24.57676
## Iteration 57 : Log-Likelihood = 24.57676##
## Call:
## lm.circular.cl(y = ..1, x = ..2, init = ..3, verbose = TRUE)
##
##
## Circular-Linear Regression
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## [1,] 5.1566 0.4184 12.325 < 2e-16 ***
## [2,] 1.1203 0.2342 4.783 8.64e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: 24.58
##
## Summary: (mu in radians)
## mu: -0.08114 ( 0.04603 ) kappa: 59.5 ( 26.5 )
## p-values are approximated using normal distributionplot(y)
lmC <- lm.circular(y=y, x=x, init=c(5,1), type='c-l', verbose=TRUE)## Iteration 1 : Log-Likelihood = 24.56042
## Iteration 2 : Log-Likelihood = 24.56847
## Iteration 3 : Log-Likelihood = 24.57257
## Iteration 4 : Log-Likelihood = 24.57465
## Iteration 5 : Log-Likelihood = 24.5757
## Iteration 6 : Log-Likelihood = 24.57622
## Iteration 7 : Log-Likelihood = 24.57649
## Iteration 8 : Log-Likelihood = 24.57662
## Iteration 9 : Log-Likelihood = 24.57669
## Iteration 10 : Log-Likelihood = 24.57672
## Iteration 11 : Log-Likelihood = 24.57674
## Iteration 12 : Log-Likelihood = 24.57675
## Iteration 13 : Log-Likelihood = 24.57675
## Iteration 14 : Log-Likelihood = 24.57675
## Iteration 15 : Log-Likelihood = 24.57675
## Iteration 16 : Log-Likelihood = 24.57675
## Iteration 17 : Log-Likelihood = 24.57676
## Iteration 18 : Log-Likelihood = 24.57676
## Iteration 19 : Log-Likelihood = 24.57676
## Iteration 20 : Log-Likelihood = 24.57676
## Iteration 21 : Log-Likelihood = 24.57676
## Iteration 22 : Log-Likelihood = 24.57676
## Iteration 23 : Log-Likelihood = 24.57676
## Iteration 24 : Log-Likelihood = 24.57676
## Iteration 25 : Log-Likelihood = 24.57676
## Iteration 26 : Log-Likelihood = 24.57676
## Iteration 27 : Log-Likelihood = 24.57676
## Iteration 28 : Log-Likelihood = 24.57676
## Iteration 29 : Log-Likelihood = 24.57676
## Iteration 30 : Log-Likelihood = 24.57676
## Iteration 31 : Log-Likelihood = 24.57676
## Iteration 32 : Log-Likelihood = 24.57676
## Iteration 33 : Log-Likelihood = 24.57676
## Iteration 34 : Log-Likelihood = 24.57676
## Iteration 35 : Log-Likelihood = 24.57676
## Iteration 36 : Log-Likelihood = 24.57676
## Iteration 37 : Log-Likelihood = 24.57676
## Iteration 38 : Log-Likelihood = 24.57676
## Iteration 39 : Log-Likelihood = 24.57676
## Iteration 40 : Log-Likelihood = 24.57676
## Iteration 41 : Log-Likelihood = 24.57676
## Iteration 42 : Log-Likelihood = 24.57676
## Iteration 43 : Log-Likelihood = 24.57676
## Iteration 44 : Log-Likelihood = 24.57676
## Iteration 45 : Log-Likelihood = 24.57676
## Iteration 46 : Log-Likelihood = 24.57676
## Iteration 47 : Log-Likelihood = 24.57676
## Iteration 48 : Log-Likelihood = 24.57676
## Iteration 49 : Log-Likelihood = 24.57676
## Iteration 50 : Log-Likelihood = 24.57676
## Iteration 51 : Log-Likelihood = 24.57676
## Iteration 52 : Log-Likelihood = 24.57676
## Iteration 53 : Log-Likelihood = 24.57676
## Iteration 54 : Log-Likelihood = 24.57676
## Iteration 55 : Log-Likelihood = 24.57676
## Iteration 56 : Log-Likelihood = 24.57676
## Iteration 57 : Log-Likelihood = 24.57676p <- circular(lmC$mu+2*atan(x%*%lmC$coefficients))
points(p, col=2, pch= "+")
The von Mises distribution, also known as the circular normal distribution, is a continuous probability distribution on the unit circle \([0, 2\pi)\). It is often used to model circular data such as angles, directions, or phases.
The probability density function (PDF) of the von Mises distribution is given by:
\[ f(\theta | \mu, \kappa) = \frac{1}{2\pi I_0(\kappa)} \exp\left(\kappa \cos(\theta - \mu)\right) \]
where: The parameter \(\kappa\) controls the shape of the distribution:The cumulative distribution function (CDF) of the von Mises distribution does not have a closed-form expression but can be expressed in terms of the modified Bessel function of the first kind.
Fitting a von Mises Distribution to Data: —————————————- —————————————
# Estimate parameters using maximum likelihood
fit <- mle.vonmises(data)
print(fit) # Prints the estimated mean direction and kappa##
## Call:
## mle.vonmises(x = data)
##
## mu: 0.6716 ( 0.08936 )
##
## kappa: 1.86 ( 0.2328 )# Plot the fitted distribution over the data
plot(data)
lines(density(data, bw = 20), col = "red", lwd = 2) 
Von Mises Density ——————-
data1 <- rvonmises(100, circular(0), 10, control.circular=list(units="degrees"))
plot(data1)
ff <- function(x) dvonmises(x, mu=circular(pi), kappa=10)
curve.circular(ff, join=TRUE, xlim=c(-2.3, 1),
main="Density of a VonMises Distribution \n mu=pi, kappa=10")
Compare with Well-Known Normal
require(graphics)
dnorm(0) == 1/sqrt(2*pi)## [1] TRUEdnorm(1) == exp(-1/2)/sqrt(2*pi)## [1] TRUEdnorm(1) == 1/sqrt(2*pi*exp(1))## [1] TRUE## Using "log = TRUE" for an extended range :
plot(function(x) dnorm(x), -5, 5,
main = "Normal density")
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
mtext("dnorm(x)", adj = 0)
mtext("dnorm(x)", col = "red", adj = 1)
Spherical Density
#install.packages("movMF") # For vMF distribution functions
#install.packages("rgl") # For 3D plotting
library(movMF)## Warning: package 'movMF' was built under R version 4.3.3library(rgl)set.seed(123)
n <- 1000 # Number of data points
p <- 3 # Dimension of the sphere (3D)
mu <- c(0, 0, 1) # Mean direction (North Pole)
kappa <- 10 # Concentration parameter
# Simulate data
data <- rmovMF(n, mu, kappa)# Create a sphere for plotting
open3d()## wgl
## 2spheres3d(0, 0, 0, radius = 1, color = "lightgray")
# Plot the data points on the sphere
points3d(data, col = "blue")
# Plot the mean direction as a red dot
points3d(mu, col = "red", size = 5)Estimate vMF Parameters:
# Estimate parameters using maximum likelihood
fit <- movMF(data, k = 1) # Fit a single vMF component
print(fit) # Prints the estimated mean direction and kappa## theta:
## [,1] [,2] [,3]
## 1 -1.439528e-05 -0.02798916 0.961476
## alpha:
## [1] 1
## L:
## [1] 141.2091Visualize Estimated Density: ——————————
# Create a grid of points on the sphere
theta <- seq(0, pi, length.out = 50)
phi <- seq(0, 2*pi, length.out = 50)
grid <- expand.grid(theta = theta, phi = phi)
# Transform to Cartesian coordinates
x <- sin(grid$theta) * cos(grid$phi)
y <- sin(grid$theta) * sin(grid$phi)
z <- cos(grid$theta)
points <- cbind(x, y, z)
# Calculate density values
density_values <- dmovMF(points, fit$theta)
# ... (previous code for simulation, fitting, and density calculation)
# Create matrices from coordinate vectors
x_mat <- matrix(x, nrow = length(theta), ncol = length(phi), byrow = TRUE)
y_mat <- matrix(y, nrow = length(theta), ncol = length(phi), byrow = TRUE)
z_mat <- matrix(z, nrow = length(theta), ncol = length(phi), byrow = TRUE)
# Plot the density as a heatmap on the sphere (corrected)
surface3d(x_mat, y_mat, z_mat,
col = heat.colors(100)[cut(density_values, 100)], alpha = 0.8)# 600 observations from two von Mises distributions
library(circular)
library(cplots)## Warning: package 'cplots' was built under R version 4.3.3x = c(rvonmises(200, circular(pi/4), 5), rvonmises(400, circular(pi), 20))
cbarplot(x)
cbarplot(x, prob=FALSE)
cbarplot(x, radius=1, nlabels=0, col="lightblue")
cbarplot(x, radius=1, col="lightblue", border="skyblue4")
2D Density
# 600 observations from two von Mises distributions
x = c(rvonmises(200, circular(pi/4), 5), rvonmises(400, circular(pi), 20))
dvm = function(x, mu=0, kappa=1) # von Mises density
exp(kappa * cos(x - mu)) * (2 * pi * besselI(kappa, 0))^(-1)
f = function(x) 1/3 * dvm(x, pi/4, 5) + 2/3 * dvm(x, pi, 20)
cdensity(f) # plot the density in an area-proportional manner
chist(x) # circular histogram
cdensity(f, add=TRUE) # superimpose the density curve
chist(x, area=FALSE) # height-proportional circular histogram
cdensity(f, area=FALSE, add=TRUE) # superimpose the density curve
chist(x, radius=0) # rose diagrams
cdensity(f, radius=0, add=TRUE)
chist(x, radius=0, area=FALSE)
cdensity(f, radius=0, area=FALSE, add=TRUE)
x = c(rvonmises(10, circular(pi/4), 5), rvonmises(20, circular(pi), 20))
cdotplot(x) 
# area-proportional dot plotcdotplot(x, area = FALSE) # height-proportional dot plot
# 900 observations from two von Mises distributions
y = c(rvonmises(300, circular(pi/4), 5), rvonmises(600, circular(pi), 20))
cdotplot(y, nbins=76, unit = 10) # area-proportional (partial) dot plot
cdotplot(y, nbins=76, unit = 10, area = FALSE) # height-proportional
# 600 observations from two von Mises distributions
x = c(rvonmises(200, circular(pi/4), 5), rvonmises(400, circular(pi), 20))
chist(x) # area-proportional circular histgram
chist(x, area = FALSE) # height-proportional circular histgram
chist(x, radius=0) # area-proportional rose diagram
chist(x, radius=0, area=FALSE) # height-proportional rose diagram
chist(x, prob=FALSE) # labels for frequency
EXERCISE
library(circular)
data("pigeons", package = "circular")x = pigeons[,2] / 180 * pi # bearing
y = pigeons[,1] # treatment
vs = split(x, factor(y, unique(y))) # list of classified value
prop = sapply(vs, length) / length(x) # proportion of each class
# Define the kde function for each class using von Mises kernels
dvm = function(x, mu=0, kappa=1) # von Mises density
exp(kappa * cos(x - mu)) * (2 * pi * besselI(kappa, 0))^(-1)
kdevm = function(x, x0, bw=0.3)
rowMeans(outer(x, x0, dvm, 0.5 / (1 - exp(-bw^2 / 2))))
fs = list(function(x) kdevm(x, x0=vs[[1]]),
function(x) kdevm(x, x0=vs[[2]]),
function(x) kdevm(x, x0=vs[[3]]))
# stacked density curves for 3 classes
cmdensity(fs) # 1:1:1
cmdensity(fs, prop)
library(circular)
data("pigeons", package = "circular")
x = pigeons[,2] / 180 * pi
y = pigeons[,1]
# stacked circular histograms
cmhist(x, y) # area-proportional
cmhist(x, y, area=FALSE) # height-proportional
library(circular)
# Simulate circular data (angles in degrees)
set.seed(123)
n <- 50
theta_uniform <- runif(n, 0, 360) # Uniformly distributed angles
theta_vonmises <- rvonmises(n, mu = circular(0), kappa = 5) # Von Mises distributed angles
# Convert to circular objects
theta_uniform_circ <- circular(theta_uniform, units = "degrees")
theta_vonmises_circ <- circular(theta_vonmises, units = "degrees")
# Visualize the data (optional)
par(mfrow=c(1,2))
plot(theta_uniform_circ, main = "Uniform Data")
plot(theta_vonmises_circ, main = "Von Mises Data")
# Rayleigh Test for Uniformity
rayleigh_uniform <- rayleigh.test(theta_uniform_circ)
rayleigh_vonmises <- rayleigh.test(theta_vonmises_circ)
print(rayleigh_uniform)##
## Rayleigh Test of Uniformity
## General Unimodal Alternative
##
## Test Statistic: 0.0296
## P-value: 0.9572print(rayleigh_vonmises)##
## Rayleigh Test of Uniformity
## General Unimodal Alternative
##
## Test Statistic: 0.9988
## P-value: 0# Watson's Test for Uniformity (and Von Mises)
watson_uniform <- watson.test(theta_uniform_circ)
watson_vonmises_uniform <- watson.test(theta_vonmises_circ) # Test for uniformity
watson_vonmises_vm <- watson.test(theta_vonmises_circ, dist = "vonmises") # Test for von Mises
print(watson_uniform)##
## Watson's Test for Circular Uniformity
##
## Test Statistic: 0.0178
## P-value > 0.10
## print(watson_vonmises_uniform)##
## Watson's Test for Circular Uniformity
##
## Test Statistic: 4.0283
## P-value < 0.01
## print(watson_vonmises_vm)##
## Watson's Test for the von Mises Distribution
##
## Test Statistic: 1.1535
## P-value < 0.01
## Circular time plots in R
Circular plots offer a powerful way to visualize and identify patterns in activities that vary throughout the day. This demonstration showcases how to use ggplot2 and circular packages in R, along with lubridate, to create insightful time-of-day visualizations."
library(lubridate)## Warning: package 'lubridate' was built under R version 4.3.3##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(ggplot2) # use at least 0.9.3 for theme_minimal()## generate random data in POSIX date-time format
set.seed(44)
N=500
events <- as.POSIXct("2011-01-01", tz="GMT") +
days(floor(365*runif(N))) +
hours(floor(24*rnorm(N))) + # using rnorm here
minutes(floor(60*runif(N))) +
seconds(floor(60*runif(N)))Next, we organize the events a little. We’re going to bin the events by hour and, to add some color, we’ll categorize events based on whether they were during the ‘workday’ (defined here as 9 thru 5pm).
# extract hour with lubridate function
hour_of_event <- hour(events)
# make a dataframe
eventdata <- data.frame(datetime = events, eventhour = hour_of_event)
# determine if event is in business hours
eventdata$Workday <- eventdata$eventhour %in% seq(9, 17)ggplot(eventdata, aes(x = eventhour, fill = Workday)) + geom_histogram(breaks = seq(0,
24), width = 2, colour = "grey") + coord_polar(start = 0) + theme_minimal() +
scale_fill_brewer() + ylab("Count") + ggtitle("Events by Time of day") +
scale_x_continuous("", limits = c(0, 24), breaks = seq(0, 24), labels = seq(0,
24))
# make our data into circular class from package circular
eventdata$eventhour <- circular(hour_of_event%%24, # convert to 24 hrs
units="hours", template="clock24")# plot a rose diagram, setting prop(ortion) argument after trial-n-error
rose.diag(eventdata$eventhour, bin = 24, col = "lightblue", main = "Events by Hour (sqrt scale)",
prop = 3)
While rose.diag() defaults to displaying the square root of values for better count visualization, we can enhance peak prominence by setting radii.scale = “linear”. Additionally, the function allows for point overlays, a feature that could be helpful when visualizing an additional variable alongside the distribution.
# redundantly add points to surface; we need to adjust parameters like
# shrink, cex, and prop
rp <- rose.diag(eventdata$eventhour, bin = 24, col = "lightblue", main = "Events by Hour (linear scale)",
prop = 12, radii.scale = "linear", shrink = 1.5, cex = 0.8)
points(eventdata$eventhour, plot.info = rp, col = "grey", stack = TRUE)
“To further analyze the distribution, we can estimate the density of the circular data and visualize it in a more familiar way. The circular.density() function is crucial here due to the cyclical nature of our angles. We’ll customize the plot for better interpretability by focusing on the shape of the curve and providing clear labels, rather than relying on numerical axis values.”
# estimate density, class is still circular not density
bw <- 10 * bw.nrd0(eventdata$eventhour) # may not be best bw: experiment
dens <- density.circular(eventdata$eventhour, bw = bw) # bw must be given# returns NULL for some reason
plot(dens, plot.type = "line", join = TRUE, main = "Probability of Event by Hour",
xlab = "Hour", ylab = "", yaxt = "n")
## NULL