Commit bc8c14be authored by Philip Rinn's avatar Philip Rinn

Clean and restore the environment in examples.r

parent 8c75a9cd
......@@ -18,10 +18,12 @@ 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)
op <- par(no.readonly = TRUE)
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 = "")
par(op)
# Esitmate drift and diffusion coefficients
ests <- Langevin1D(Ux, bins, steps)
......@@ -29,6 +31,7 @@ attach(ests)
# Plot drift and diffusion coefficients with errors and theoretical expecations
# for comparison
op <- par(no.readonly = TRUE)
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.]",
......@@ -39,11 +42,12 @@ plotCI(mean_bin, D2, uiw = eD2, xlab = "x [a.u.]",
ylab = expression(paste("Diffusion coefficient ", D^(2), "(x) [",1/s^2,"]")),
pch = 20)
lines(mean_bin, mean_bin^2 + 1, col = 'red', lty = 2)
par(op)
# 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))
detach(ests)
# Generate a time series from the reconstructed coefficients
set.seed(5048);
......@@ -53,6 +57,7 @@ rec_x <- timeseries1D(N = 1e7, d10 = estD1[1], d11 = estD1[2], d12 = estD1[3],
# Compute and plot the increment PDFs for the original and the reconstructed
# time series for four different time lags
op <- par(no.readonly = TRUE)
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),
......@@ -69,6 +74,7 @@ for(i in 1:4) {
lines(den, col = i)
lines(rec_den, lty = 2, col = i)
}
par(op)
##########################################
......@@ -105,12 +111,14 @@ 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
op <- par(no.readonly = TRUE)
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')
par(op)
# Estimate drift and diffusion coefficients
ests <- Langevin2D(Ux, bins, steps)
......@@ -121,6 +129,7 @@ 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
op <- par(no.readonly = TRUE)
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))
......@@ -173,3 +182,5 @@ 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)
par(op)
detach(ests)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment