R: Example scripts and user provided functions for the Gaussian...
examples
R Documentation
Example scripts and user provided functions for the Gaussian Diffusion Process
Description
For a generic Gaussian Diffusion process X(t) the drift can be written as a(t)*X(t)+b(t) and the infinitesimal variance as cc(t)^2. User should provide the functional form for a(t),b(t) and cc(t) in the main script. The threshold to be crossed has to be also provided through the function S(t) and its derivative Sp(t) in the same script. Examples of such a script are given for a Wiener process with or without drift through different thresholds (scripts Wiener, Wiener1, WienerDrift) for an Ornstein Uhlenbeck process also in presence of an additional current (scripts OrnUhl, OrnUhlCurrent) and for a more complex process with time-varying coefficients (Logistic).
Usage
Wiener.R
Author(s)
A. Buonocore, M.F. Carfora
Examples
##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
# Wiener process through a periodic boundary: the process is built, its FPT pdf
# evaluated by numerical quadrature and a train of crossing times is simulated.
# Results are shown and saved to a file
library(GaDiFPT)
cat('#################################################################### \n')
cat('##### First Passage Time Simulation ##### \n')
cat('##### for the Wiener process ##### \n')
cat('##### through a periodic boundary ##### \n')
cat('#################################################################### \n \n \n')
Wiener1 <- diffusion(c("mu","sigma2"))
mu <- 0.0
sigma2 <- 1.0
S0 <- 10
S1 <- 3
Sfr <- 0.005
Stype <- "periodic"
t0 <- 0.0
x0 <- 0.0
Tfin <- 1000
deltat <- 0.5
N <- floor((Tfin - t0)/deltat)
M <- 3000
quadflag <- 1
RStudioflag <- TRUE
param <- inputlist(mu,sigma2,Stype,t0,x0,Tfin,deltat,M,quadflag,RStudioflag)
fileout <- "results_Wiener1.out"
aaa <- function(t) {
aaa <- 0.0 + 0.0*t
}
bbb <- function(t) {
bbb <- mu + 0.0*t
}
SSS <- function(t) {
SSS <- S0+S1*cos(Sfr*t)
}
SSSp <- function(t) {
SSSp <- -S1*Sfr*sin(Sfr*t)
}
mp <- numeric(N+1)
up <- numeric(N+1)
vp <- numeric(N+1)
app <- numeric(N)
tempi <- seq(t0, by=deltat, length=N+1)
dum <- vectorsetup(param)
mp <- dum[,1]
up <- dum[,2]
vp <- dum[,3]
splot <- S(tempi)
mp1 <- mp - sqrt(2*sigma2)
mp2 <- mp + sqrt(2*sigma2)
matplot(tempi, cbind(mp,mp1,mp2,splot),type="l",lty=c(1,2,2,1),lwd=1,
main="mean of the process vs. threshold",xlab="time(ms)",ylab="")
legend("bottomright",c("mean","threshold"),
lty=c(1,1),col=c("black","blue"))
Nmax <- which.min(abs(mp[2:(N+1)]-mp[1:N]))
N1 <- N
if (quadflag == 0) N1 <- max(c(Nmax,N/4))
N1p1 <- N1+1
answer <- FPTdensity_byint(param,N1)
plot(answer)
spikes <- FPTsimul(answer,M)
histplot(spikes,answer)
res_summary(answer,M,fileout)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(GaDiFPT)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GaDiFPT/examples.Rd_%03d_medium.png", width=480, height=480)
> ### Name: examples
> ### Title: Example scripts and user provided functions for the Gaussian
> ### Diffusion Process
> ### Aliases: examples Wiener Wiener1 WienerDrift OrnUhl OrnUhlCurrent
> ### Logistic
>
> ### ** Examples
>
> ##---- Should be DIRECTLY executable !! ----
> ##-- ==> Define data, use random,
> ##-- or do help(data=index) for the standard data sets.
>
> # Wiener process through a periodic boundary: the process is built, its FPT pdf
> # evaluated by numerical quadrature and a train of crossing times is simulated.
> # Results are shown and saved to a file
>
>
> library(GaDiFPT)
>
> cat('#################################################################### \n')
####################################################################
> cat('##### First Passage Time Simulation ##### \n')
##### First Passage Time Simulation #####
> cat('##### for the Wiener process ##### \n')
##### for the Wiener process #####
> cat('##### through a periodic boundary ##### \n')
##### through a periodic boundary #####
> cat('#################################################################### \n \n \n')
####################################################################
>
> Wiener1 <- diffusion(c("mu","sigma2"))
>
> mu <- 0.0
> sigma2 <- 1.0
>
> S0 <- 10
> S1 <- 3
> Sfr <- 0.005
> Stype <- "periodic"
>
> t0 <- 0.0
> x0 <- 0.0
> Tfin <- 1000
> deltat <- 0.5
> N <- floor((Tfin - t0)/deltat)
> M <- 3000
>
> quadflag <- 1
> RStudioflag <- TRUE
>
> param <- inputlist(mu,sigma2,Stype,t0,x0,Tfin,deltat,M,quadflag,RStudioflag)
>
> fileout <- "results_Wiener1.out"
>
> aaa <- function(t) {
+ aaa <- 0.0 + 0.0*t
+ }
>
> bbb <- function(t) {
+ bbb <- mu + 0.0*t
+ }
>
> SSS <- function(t) {
+ SSS <- S0+S1*cos(Sfr*t)
+ }
>
> SSSp <- function(t) {
+ SSSp <- -S1*Sfr*sin(Sfr*t)
+ }
>
> mp <- numeric(N+1)
> up <- numeric(N+1)
> vp <- numeric(N+1)
> app <- numeric(N)
>
> tempi <- seq(t0, by=deltat, length=N+1)
>
> dum <- vectorsetup(param)
> mp <- dum[,1]
> up <- dum[,2]
> vp <- dum[,3]
>
> splot <- S(tempi)
> mp1 <- mp - sqrt(2*sigma2)
> mp2 <- mp + sqrt(2*sigma2)
> matplot(tempi, cbind(mp,mp1,mp2,splot),type="l",lty=c(1,2,2,1),lwd=1,
+ main="mean of the process vs. threshold",xlab="time(ms)",ylab="")
> legend("bottomright",c("mean","threshold"),
+ lty=c(1,1),col=c("black","blue"))
>
>
> Nmax <- which.min(abs(mp[2:(N+1)]-mp[1:N]))
>
> N1 <- N
> if (quadflag == 0) N1 <- max(c(Nmax,N/4))
>
> N1p1 <- N1+1
>
> answer <- FPTdensity_byint(param,N1)
> plot(answer)
>
> spikes <- FPTsimul(answer,M)
> histplot(spikes,answer)
>
> res_summary(answer,M,fileout)
####################################################################
##### First Passage Time Simulation #####
##### for Gaussian Diffusions: Results #####
####################################################################
PROBABILITY MASS REACHED = 0.7377655
Mean time from the density (ms) = 193.2596099
Standard deviation (ms) = 220.1627624
Median time from the distribution (ms) = 304.5
Mean time for the spikes (ms) = 1118.761194
Standard deviation (ms) = 2002.871734
Median time for the spikes (ms) = 294.1322275
>
>
>
>
>
> dev.off()
null device
1
>