Skip to contents

Simulate a blurred single or multichannel spike train

Source: R/simulate_spike_train.R
simulate_spike_train.Rd

Generate a sparse spike train blurred by a known point spread function (PSF), with Poisson, Gaussian or negative binomial noise. Columns of the design are shifted copies of the PSF. This models signal deconvolution, including GC/MS data with sparse mass spectra shared over time.

Usage

simulate_spike_train(
  n = 100,
  p = n,
  k = 10,
  mean_beta = 10000,
  sd_logbeta = 1,
  family = "poisson",
  sd_noise = 10,
  psf = function(x, u, w = 2.5) exp((x - u)^2/(-2 * w^2)),
  thresh = 1e-05,
  sparse = TRUE,
  multichannel = FALSE,
  seed = 123,
  Plot = FALSE,
  intercept = FALSE,
  theta = NULL,
  m = NULL,
  spectra = NULL,
  channel.sparsity = 0.7,
  clip = FALSE
)

Arguments

n, p

Numbers of time points and predictor locations; p>n allows superresolution. The design excludes any intercept column.

k

Number of peaks, capped at min(n,p).

mean_beta

Median total peak height (the log-normal scale parameter).

sd_logbeta

Standard deviation of log peak heights. Heights follow exp(rnorm(k, log(mean_beta), sd_logbeta)); their arithmetic expectation is mean_beta*exp(sd_logbeta^2/2).

family

"poisson", "gaussian", "nb" or "negative.binomial", or a corresponding family object with an identity link. NB family objects are supplied by mgcv::nb(link="identity"). Identity-link signal_families objects also support fixed-parameter NB1, common-scale Gamma and censored Gamma noise.

sd_noise

Gaussian noise standard deviation, one value or one per channel.

psf

Function of time x and peak location u. The default Gaussian has height one and width 2.5. Values must be finite and nonnegative.

thresh

PSF entries smaller than this nonnegative cutoff are zeroed.

sparse

Return a sparse design/coefficient matrix? Defaults to TRUE. The sparse design is constructed directly without a dense n by p copy.

multichannel

Generate multiple spectral channels? With m=NULL, use the package's simulated_spectra library, retaining informative channels.

seed

Random seed; a specified seed preserves the caller's RNG state. NULL uses and advances the current stream.

Plot

Plot the generated signal or multichannel heatmap?

intercept

FALSE/NULL omits an intercept; TRUE adds a baseline of one; a scalar or channel vector supplies its value, including zero. The corresponding row is included in beta_true, but not in X.

theta

Positive NB shape, one value or one per channel. Defaults to 10 or the supplied family's fixed shape. NB variance is mu+mu^2/theta.

m

Explicit number of channels. Values greater than one enable multichannel generation using random sparse spectra, unless spectra is provided. NULL selects the library when multichannel=TRUE.

spectra

Optional nonnegative matrix of template spectra, with one spectrum per row. Rows are normalized to sum to one, then k rows are sampled, with replacement only if the library contains fewer than k rows. Channels are retained, including channels with zero generated signal.

channel.sparsity

Probability of a zero entry in random spectra (m supplied); defaults to 0.7. Each peak retains at least one channel.

clip

Clip negative Gaussian observations to zero? Defaults to FALSE for untruncated Gaussian noise. Clipping gives a left-censored response. Numeric simul() setups 1/3 enable clipping by default.

Value

A list with X, beta_true, y, y_true, eta_true, family, weights, intercept, intercept_true, dispersion, theta, spectra, time grid x, and predictor-location grid x_beta. Means satisfy y_true = X %*% slopes + intercept_true. beta_true includes the intercept row when requested. Use intercept=sim$intercept when fitting sim$X. With clip=TRUE, y_true describes the latent uncensored signal.

See also

simul(), L0glm.fit(), sparseglm()