sinewavePDF: PDF of a truncated sinusoidal curve

In ADMUR: Ancient Demographic Modelling Using Radiocarbon

Description

Probability density function for a truncated sinusoidal curve. \loadmathjax

Usage

sinewavePDF(x, min, max, f, p, r)

Arguments

x	Numeric vector of years
min, max	Lower and upper \mjseqnx limits of the distribution
f	Numeric frequency (cycles per unit \mjseqnx).
p	Numeric between \mjseqn0 and \mjseqn2\pi, giving the cycle position (in radians) at \mjseqnx = 0.
r	Numeric between 0 and 1, determining how flat the distribution is.

Details

The usual function to describe a sine wave is \mjseqnf(x) = A\sin(2\pi f x + p), where \mjseqnA is the amplitude, \mjseqnf is the frequency (cycles per year), and \mjseqnp is the cycle position (in radians) at \mjseqnx = 0, and therefore oscillates above and below the x-axis.

However, a sinusoidal PDF must by definition always be non-negative, which can conceptually be considered as a sine wave stacked on top of a uniform distribution with a height \mjseqnA + k, where \mjseqnk >= 0. Since the PDF is \mjseqnf(x) divided by the area below the curve, A and k simplify to a single parameter \mjseqnr that determines the relative proportions of the uniform and sinusoidal components, such that:

when \mjseqnr = 0 the amplitude of the sine wave component is zero, and the overall PDF is just a uniform distribution between min and max.

when \mjseqnr = 1 the uniform component is zero, and the minima of the sine wave touches zero. This does not necessarily mean the PDF minimum equals zero, since a minimum point of the sine wave may not occur with PDF domain (truncated between min and max).

Therefore the formula for the PDF is:

\mjsdeqn\frac

1 + \sin(2\pi f x + p) - \ln(r)(x_max - x_min)(1 - \ln(r)) + (\frac12\pi f)[\cos(2\pi f x_min - p) - \cos(2\pi f x_max - p)] where \mjseqnx = years, and \mjseqnx_min and \mjseqnx_max determine the truncated date range.

Examples

	# A sinewave with a period of 700 years
	x <- seq(1500,4500, length.out=1000)
	y <- sinewavePDF(x, min=2000, max=4000, f=1/700, p=0, r=0.2)
	plot(x,y,type='l')
	

ADMUR documentation built on March 24, 2021, 1:08 a.m.