examples.r 7.13 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
library("Langevin")
# Load library used for plotting
library("plotrix")

# Set the number of bins to estimate drift and diffusion on
bins <- 40
# Set the time lags for computing the conditional moments
steps <- c(1:3)
# Set the sampling frequency
sf <- 1000

##########################################
# One dimensional example
##########################################

# Generate the time series
set.seed(4711)
Ux <- timeseries1D(N = 1e7, d11 = 1, d13 = -1, d22 = 1, d20 = 1, sf = sf)

# Plot time series and probability density function (PDF)
21
op <- par(no.readonly = TRUE)
22 23 24 25
par(mar = c(4.2, 5.4, 0.5, 0.5))
layout(matrix(c(1, 1, 2), 1, 3))
plot((1:1e5)/sf, Ux[1:1e5], xlab = "t [s]", ylab = "x [a.u.]", t = 'l')
plot(density(Ux), xlab = "x [a.u.]", ylab = "density", t = 'l',main = "")
26
par(op)
27 28 29 30 31 32 33

# Esitmate drift and diffusion coefficients
ests <- Langevin1D(Ux, bins, steps)
attach(ests)

# Plot drift and diffusion coefficients with errors and theoretical expecations
# for comparison
34
op <- par(no.readonly = TRUE)
35 36 37 38 39 40 41
par(mar = c(4.2, 5.4, 0.5, 0.5))
par(mfrow = c(1, 2))
plotCI(mean_bin, D1, uiw = eD1, xlab = "x [a.u.]",
       ylab = expression(paste("Drift coefficient ", D^(1), "(x) [",1/s,"]")),
       pch = 20)
lines(mean_bin, mean_bin - mean_bin^3, col = 'red', lty = 2)
plotCI(mean_bin, D2, uiw = eD2, xlab = "x [a.u.]",
42
       ylab = expression(paste("Diffusion coefficient ", D^(2), "(x) [",1/s,"]")),
43 44
       pch = 20)
lines(mean_bin, mean_bin^2 + 1, col = 'red', lty = 2)
45
par(op)
46 47 48 49

# Fit polynomials to the estimated drift and diffusion coefficients
estD1 <- coef(lm(D1 ~ mean_bin + I(mean_bin^2) + I(mean_bin^3), weights = 1/eD1))
estD2 <- coef(lm(D2 ~ mean_bin + I(mean_bin^2), weights = 1/eD2))
50
detach(ests)
51 52 53 54 55 56 57 58 59

# Generate a time series from the reconstructed coefficients
set.seed(5048);
rec_x <- timeseries1D(N = 1e7, d10 = estD1[1], d11 = estD1[2], d12 = estD1[3],
                      d13 = estD1[4], d20 = estD2[1], d21 = estD2[2],
                      d22 = estD2[3], sf = sf)

# Compute and plot the increment PDFs for the original and the reconstructed
# time series for four different time lags
60
op <- par(no.readonly = TRUE)
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
par(mar = c(4.2, 5.4, 0.5, 0.5))
par(mfrow = c(1, 1));
plot(1, 1, log = "y", type = "n", xlim = c(-11, 12), ylim = c(1e-8, 5),
     xlab = expression(Delta[tau]/sigma[Delta[tau]]), ylab = "density")

tau <- c(1, 10, 100, 1000)
for(i in 1:4) {
    delta <- diff(Ux, lag = tau[i])
    rec_delta <- diff(rec_x, lag = tau[i])
    den <- density(delta)
    den$x <- den$x/sd(delta, na.rm = TRUE)
    rec_den <- density(rec_delta)
    rec_den$x <- rec_den$x/sd(rec_delta, na.rm = TRUE)
    lines(den, col = i)
    lines(rec_den, lty = 2, col = i)
}
77
par(op)
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113


##########################################
# Two dimensional example
##########################################

# Construct the drift matrices
D1_1 <- matrix(0, nrow = 4,ncol = 4)
D1_1[1, 2]  =  1
D1_2 <- matrix(0, nrow = 4,ncol = 4)
D1_2[2, 1]  =  0.02
D1_2[4, 1]  =  -1
D1_2[1, 2]  =  0.03
D1_2[3, 2]  =  -1

# Construct the diffusion matrices
g_11 <- matrix(0, nrow = 3, ncol = 3)
g_12 <- matrix(0, nrow = 3, ncol = 3)
g_21 <- matrix(0, nrow = 3, ncol = 3)
g_22 <- matrix(0, nrow = 3, ncol = 3)

# Generate the time series without noise (only for plotting)
set.seed(4711)
Ux1 <- timeseries2D(N = 1e6, 0.145, 0.0002, D1_1, D1_2, g_11, g_12, g_21,
                    g_22, sf = sf)

# Add the noise
g_11[1, 1] <- 0.0025
g_22[1, 1] <- 0.0025

# Generate time series with added noise
set.seed(4711)
Ux <- timeseries2D(N = 1e8, 0.145, 0.0002, D1_1, D1_2, g_11, g_12, g_21,
                   g_22, sf = sf)

# Plot trajectories of the time series without and with noise
114
op <- par(no.readonly = TRUE)
115 116 117 118 119 120
par(mar = c(4.2, 5.4, 0.5, 0.5))
par(mfrow = c(1, 2))
plot(Ux1[1:1e6, 1], Ux1[1:1e6, 2], xlab = expression(paste(X[1]," [a.u.]")),
     ylab = expression(paste(X[2]," [a.u.]")), t = 'l')
plot(Ux[1:1e5, 1], Ux[1:1e5, 2], xlab = expression(paste(X[1]," [a.u.]")),
     ylab = expression(paste(X[2]," [a.u.]")), t = 'l')
121
par(op)
122 123 124 125 126 127 128 129 130 131

# Estimate drift and diffusion coefficients
ests <- Langevin2D(Ux, bins, steps)
attach(ests)
# Compute the mids of the bins for plotting
bin_mids <- matrix(NA, nrow = 40, ncol = 2)
bin_mids[,1] <- U[1:40, 1] + diff(U[,1])/2
bin_mids[,2] <- U[1:40, 2] + diff(U[,2])/2

# Plot drift and diffusion coefficients
132
op <- par(no.readonly = TRUE)
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
par(mar = c(4.2, 5.4, 0.5, 0.5))
layout(matrix(c(1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5), 2, 6, byrow = TRUE))

pmat <- perspx(x = bin_mids[,1], y = bin_mids[,2], z = D1[,,1], xlab = "", ylab = "",
               zlab = "", ticktype = "detailed", theta = 30, main = "")
attr(pmat, "ranges") <- list(x = range(bin_mids[,1]), y = range(bin_mids[,2]),
                             z = range(D1[,,1], na.rm = TRUE))
mtext3d("X", pmat, expression(X[1]))
mpos<-plotrix:::get_axispos3d("Y", pmat, at = 0, dist = .3)
mpos[1,] <- -mpos[1,]
ptext3d(mpos[1,], mpos[2,], mpos[3,], pmat = pmat, xpd = NA, texts = expression(X[2]))
mtext3d("Z", pmat, expression(D[x]^(1)), at = 0)

pmat <- perspx(x = bin_mids[,1], y = bin_mids[,2], z = D1[,,2], xlab = "", ylab = "",
               zlab = "", ticktype = "detailed", theta = 30, main = "")
attr(pmat, "ranges") <- list(x = range(bin_mids[,1]), y = range(bin_mids[,2]),
                             z = range(D1[,,2], na.rm = TRUE))
mtext3d("X", pmat, expression(X[1]))
mpos<-plotrix:::get_axispos3d("Y", pmat, at = 0, dist = .3)
mpos[1,] <- -mpos[1,]
ptext3d(mpos[1,], mpos[2,], mpos[3,], pmat = pmat, xpd = NA, texts = expression(X[2]));
mtext3d("Z", pmat, expression(D[y]^(1)), at = 0)

pmat <- perspx(x = bin_mids[,1], y = bin_mids[,2], z = D2[,,1], xlab = "", ylab = "",
               zlab = "", ticktype = "detailed", theta = 30, main = "")
attr(pmat, "ranges") <- list(x = range(bin_mids[,1]), y = range(bin_mids[,2]),
                             z = range(D2[,,1], na.rm = TRUE))
mtext3d("X", pmat, expression(X[1]))
mpos<-plotrix:::get_axispos3d("Y", pmat, at = 0, dist = .3)
mpos[1,] <- -mpos[1,]
ptext3d(mpos[1,], mpos[2,], mpos[3,], pmat = pmat, xpd = NA, texts = expression(X[2]));
mtext3d("Z", pmat, expression(D[xx]^(2)), at = 0.0025)

pmat <- perspx(x = bin_mids[,1], y = bin_mids[,2], z = D2[,,3], xlab = "", ylab = "",
               zlab = "", ticktype = "detailed", theta = 30, main = "")
attr(pmat, "ranges") <- list(x = range(bin_mids[,1]), y = range(bin_mids[,2]),
                             z = range(D2[,,3], na.rm = TRUE))
mtext3d("X", pmat, expression(X[1]))
mpos<-plotrix:::get_axispos3d("Y", pmat, at = 0, dist = .3)
mpos[1,] <- -mpos[1,]
ptext3d(mpos[1,], mpos[2,], mpos[3,], pmat = pmat, xpd = NA, texts = expression(X[2]));
mtext3d("Z", pmat, expression(D[yy]^(2)), at = 0.00225)

pmat <- perspx(x = bin_mids[,1], y = bin_mids[,2], z = D2[,,2], xlab = "", ylab = "",
               zlab = "", ticktype = "detailed", theta = 30, main = "")
attr(pmat, "ranges") <- list(x = range(bin_mids[,1]), y = range(bin_mids[,2]),
                             z = range(D2[,,2], na.rm = TRUE))
mtext3d("X", pmat, expression(X[1]))
mpos<-plotrix:::get_axispos3d("Y", pmat, at = 0, dist = .3)
mpos[1,] <- -mpos[1,]
ptext3d(mpos[1,], mpos[2,], mpos[3,], pmat = pmat, xpd = NA, texts = expression(X[2]))
mtext3d("Z", pmat, expression(D[xy]^(2)), at = -2.5e-4)
185 186
par(op)
detach(ests)