skedastic (version 1.0.0)

dpeak: Probability mass function of number of peaks in an i.i.d. random sequence

Description

This function computes \(P(n,k)\) as defined by Goldfeld65;textualskedastic, i.e. the probability that a sequence of \(n\) independent and identically distributed random variables contains exactly \(k\) peaks, with peaks also as defined by Goldfeld65;textualskedastic. The function is used in ppeak to compute \(p\)-values for the Goldfeld-Quandt nonparametric test for heteroskedasticity in a linear model.

Usage

dpeak(k, n, usedata = FALSE)

Arguments

k

An integer or a sequence of integers strictly incrementing by 1, with all values between 0 and n - 1 inclusive. Represents the number of peaks in the sequence.

n

A positive integer representing the number of observations in the sequence.

usedata

A logical. Should probability mass function values be read from dpeakdat rather than computing them? This option will save significantly on computation time if \(n < 170\) but is currently only available for \(n \leq 1000\).

Value

A double between 0 and 1 representing the probability of exactly k peaks occurring in a sequence of \(n\) independent and identically distributed continuous random variables. The double has a names attribute with the values corresponding to the probabilities. Computation time is very slow for \(n > 170\) (if usedata is FALSE) and for \(n > 1000\) regardless of usedata value.

References

See Also

ppeak, goldfeld_quandt

Examples

# NOT RUN {
dpeak(0:9, 10)
plot(0:9, dpeak(0:9, 10), type = "p", pch = 20, xlab = "Number of Peaks",
         ylab = "Probability")

# This would be extremely slow if usedata were set to FALSE:
prob <- dpeak(0:199, 200, usedata = TRUE)
plot(0:199, prob, type = "l", xlab = "Number of Peaks", ylab = "Probability")

# `dpeakdat` is a dataset containing probabilities generated from `dpeak`
utils::data(dpeakdat)
expval <- unlist(lapply(dpeakdat,
                 function(p) sum(p * 0:(length(p) - 1))))
plot(1:1000, expval[1:1000], type = "l", xlab = parse(text = "n"),
     ylab = "Expected Number of Peaks")
# }