contMC: Continuous time Markov chain

In stepR: Multiscale Change-Point Inference

Continuous time Markov chain

Description

Simulate a continuous time Markov chain.

Deprecation warning: This function is mainly used for patchlamp recordings and may be transferred to a specialised package.

Usage

contMC(n, values, rates, start = 1, sampling = 1, family = c("gauss", "gaussKern"),
  param = NULL)

Arguments

n	number of data points to simulate
values	a numeric vector specifying signal amplitudes for different states
rates	a square matrix matching the dimension of values each with rates[i,j] specifying the transition rate from state i to state j; the diagonal entries are ignored
start	the state in which the Markov chain is started
sampling	the sampling rate
family	whether Gaussian white ("gauss") or coloured ("gaussKern"), i.e. filtered, noise should be added; cf. family
param	for family="gauss", a single non-negative numeric specifying the standard deviation of the noise; for family="gaussKern", param must be a list with entry df giving the dfilter object used for filtering, an integer entry over which specifies the oversampling factor of the filter, i.e. param$df has to be created for a sampling rate of sampling times over, and an additional non-negative numeric entry sd specifying the noise's standard deviation after filtering; cf. family

Value

A list with components

cont	an object of class stepblock containing the simulated true values in continuous time, with an additional column state specifying the corresponding state
discr	an object of class stepblock containing the simulated true values reduced to discrete time, i.e. containing only the observable blocks
data	a data.frame with columns x and y containing the times and values of the simulated observations, respectively

Note

This follows the description for simulating ion channels given by VanDongen (1996).

References

VanDongen, A. M. J. (1996) A new algorithm for idealizing single ion channel data containing multiple unknown conductance levels. Biophysical Journal 70(3), 1303–1315.

See Also

stepblock, jsmurf, stepbound, steppath, family, dfilter

Examples

# Simulate filtered ion channel recording with two states
set.seed(9)
# sampling rate 10 kHz
sampling <- 1e4
# tenfold oversampling
over <- 10
# 1 kHz 4-pole Bessel-filter, adjusted for oversampling
cutoff <- 1e3
df <- dfilter("bessel", list(pole=4, cutoff=cutoff / sampling / over))
# two states, leaving state 1 at 1 Hz, state 2 at 10 Hz
rates <- rbind(c(0, 1e0), c(1e1, 0))
# simulate 5 s, level 0 corresponds to state 1, level 1 to state 2
# noise level is 0.1 after filtering
sim <- contMC(5 * sampling, 0:1, rates, sampling=sampling, family="gaussKern",
  param = list(df=df, over=over, sd=0.1))
sim$cont
plot(sim$data, pch = ".")
lines(sim$discr, col = "red")
# noise level after filtering, estimated from first block
sd(sim$data$y[1:sim$discr$rightIndex[1]])
# show autocovariance in first block
acf(ts(sim$data$y[1:sim$discr$rightIndex[1]], freq=sampling), type = "cov")
# power spectrum in first block
s <- spec.pgram(ts(sim$data$y[1:sim$discr$rightIndex[1]], freq=sampling), spans=c(200,90))
# cutoff frequency is where power spectrum is halved
abline(v=cutoff, h=s$spec[1] / 2, lty = 2)

stepR documentation built on Nov. 14, 2023, 1:09 a.m.