smash | R Documentation |
This is a wrapper function for
smash.gaus
or smash.poiss
as
appropriate. For details see smash.gaus
and
smash.poiss
.
smash(x, model = NULL, ...)
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. |
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.
See smash.gaus
or smash.poiss
for details.
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.