smashr (version 1.3-12)

smash: Estimate the underlying mean or intensity function from Gaussian or Poisson data, respectively.

Description

This is a wrapper function for smash.gaus or smash.poiss as appropriate. For details see smash.gaus and smash.poiss.

Usage

smash(x, model = NULL, ...)

Value

See smash.gaus or smash.poiss

for details.

Arguments

x

A vector of observations (assumed equally spaced).

model

Specifies the model (Gaussian or Poisson). Can be NULL, in which case the Poisson model is assumed if x consists of integers, and the Gaussian model is assumed otherwise. One of 'gaus' or 'poiss' can also be specified to fit a specific model.

...

Additional arguments passed to smash.gaus or smash.poiss.

Details

Performs nonparametric regression on univariate Poisson or Gaussian data using wavelets. Unlike many wavelet methods the data do not need to have length that is a power of 2 or be cyclic -- the functions internally deal with these issues. For the Poisson case, the data are assumed to be i.i.d. from an underlying inhomogeneous mean function that is "smooth". Similarly for the Gaussian case, the data are assumed to be independent with an underlying smooth mean function. In the Gaussian case, the variances are allowed vary, but are similarly "spatially structured" as with the mean function. The functions smash.gaus and smash.pois perform smoothing for Gaussian and Poisson data respectively.

Examples

# First Gaussian example
# ----------------------
# In this first example, the length of the signal is a power of 2.
#
# Create the baseline mean function. (The "spikes" function is used
# as an example here.)
set.seed(2)
n <- 2^9
t <- 1:n/n
spike.f <- function (x) (0.75 * exp(-500 * (x - 0.23)^2) +
  1.5 * exp(-2000 * (x - 0.33)^2) + 3 * exp(-8000 * (x - 0.47)^2) +
  2.25 * exp(-16000 * (x - 0.69)^2) + 0.5 * exp(-32000 * (x - 0.83)^2))
mu.s <- spike.f(t)

# Scale the signal to be between 0.2 and 0.8
mu.t <- (1 + mu.s)/5
plot(mu.t,type = "l")

# Create the baseline variance function. (The function V2 from Cai &
# Wang (2008) is used here.)
var.fn <- (1e-04 + 4 * (exp(-550 * (t - 0.2)^2) +
                        exp(-200 * (t - 0.5)^2) +
                        exp(-950 * (t - 0.8)^2)))/1.35
plot(var.fn,type = "l")

# Set the signal-to-noise ratio.
rsnr    <- sqrt(5)
sigma.t <- sqrt(var.fn)/mean(sqrt(var.fn)) * sd(mu.t)/rsnr^2

# Simulate an example dataset.
X.s <- rnorm(n,mu.t,sigma.t)

# Run smash (Gaussian version is run since observations are not
# counts).
mu.est <- smash(X.s)

# Plot the true mean function as well as the estimated one.
plot(mu.t,type = "l")
lines(mu.est,col = 2)

# First Poisson example
# ---------------------
# Scale the signal to be non-zero and to have a low average intensity.
mu.t <- 0.01 + mu.s

# Simulate an example dataset.
X.s <- rpois(n,mu.t)

# Run smash (the Poisson version is run since observations are counts).
mu.est <- smash(X.s)

# Plot the true mean function as well as the estimated one.
plot(mu.t,type = "l")
lines(mu.est,col = 2)

# Second Gaussian example
# -----------------------
# In this second example, we demonstrate that smash also works even
# when the signal length (n) is not a power of 2.
n    <- 1000
t    <- 1:n/n
mu.s <- spike.f(t)

# Scale the signal to be between 0.2 and 0.8.
mu.t <- (1 + mu.s)/5

# Create the baseline variance function.
var.fn  <- (1e-04 + 4 * (exp(-550 * (t - 0.2)^2) +
                         exp(-200 * (t - 0.5)^2) +
                         exp(-950 * (t - 0.8)^2)))/1.35
sigma.t <- sqrt(var.fn)/mean(sqrt(var.fn)) * sd(mu.t)/rsnr^2

# Simulate an example dataset.
X.s <- rnorm(n,mu.t,sigma.t)

# Run smash.
mu.est <- smash(X.s)

# Plot the true mean function as well as the estimated one.
plot(mu.t,type = "l")
lines(mu.est,col = 2)

# Second Poisson example
# ----------------------
# The Poisson version of smash also works with signals that are not
# exactly of length 2^J for some integer J.
#
# Scale the signal to be non-zero and to have a low average intensity.
mu.t <- 0.01 + mu.s

# Simulate an example dataset
X.s <- rpois(n,mu.t)

# Run smash (Poisson version is run since observations are counts).
mu.est <- smash(X.s)

# Plot the true mean function as well as the estimated one.
plot(mu.t,type = "l")
lines(mu.est,col = 2)