contMC: 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

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

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.

Examples

Run this code

# NOT RUN {
# 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)
# }

Run the code above in your browser using DataLab

Data engineering and BI courses are free!