Commit 8c75a9cd authored by Philip Rinn's avatar Philip Rinn

Add source code of examples and make GitLab happy with a LICENSE file

parent 7b8b9dbd
^.*\.Rproj$
^\.Rproj\.user$
R/RcppExports.R
LICENSE
examples.r
GPL-2
\ No newline at end of file
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)
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 = "")
# Esitmate drift and diffusion coefficients
ests <- Langevin1D(Ux, bins, steps)
attach(ests)
# Plot drift and diffusion coefficients with errors and theoretical expecations
# for comparison
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.]",
ylab = expression(paste("Diffusion coefficient ", D^(2), "(x) [",1/s^2,"]")),
pch = 20)
lines(mean_bin, mean_bin^2 + 1, col = 'red', lty = 2)
# 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))
# 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
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)
}
##########################################
# 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
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')
# 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
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)
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