Phase Function
About Online computation C++ Class Mathematica Package Matlab function R function
R function
Listing
 1 phasefunc <- function(A, B, tau) {
 2   # Computes phase function at the specified lag (tau) and returns it.
 3   #
 4   # Usage:
 5   #       phasefunc(A, B, tau)
 6   #
 7   # Args:
 8   #       A: a vector holding spike/discrete time events,
 9   #       B: another vector (A itself for finding auto phase function),
10   #       tau: lag between the vector at which the phase function is sought.
11   #
12   # Returns:
13   #       Phase function between A and B vectors computed at lag tau.
14   #
15   # Reference:
16   #       A Phase Function to Quantify Serial Dependence between Discrete Samples.
17   #       Ramana Dodla and Charles J. Wilson.
18   #       Biophys J. 2010 February 17; 98(4): L5-L7.
19   #       http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2820633/
20   #
21   # Example: The following code plots a curve in Fig. 2b of the above reference.
22   #
23   #    A <- c(0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024)
24   #    tau <- seq(0, 0.4, 0.001)
25   #    psi <- as.vector(numeric(0))
26   #    for(i in seq(1, length(tau), 1)) {
28   #      psi <- c(psi, phasefunc(A, A, tau[i]))
29   #    }
30   #    plot(tau, psi)
31   #
32   # Author: Ramana Dodla, 2011
33   #
34   tau <- as.numeric(tau)
35   if(length(tau) != 1) stop("length of tau (i.e. lag, the third argument) must be 1")
36   A <- as.numeric(as.vector(A))
37   if(length(A) < 3) stop("length of A must be > 2 ")
38   B <- as.numeric(as.vector(B))
39   if(length(B) < 3) stop("length of B must be > 2 ")
40   B <- B + tau
41   # --------------------------------------
42   # Compute (gamma_i, delta_i) phase pairs
43   # --------------------------------------
44   i=2
45   j=2
46   pairs <- matrix(numeric(0), 0, 2)
47   Tval <- max(A[1],B[1])
48   while (i <= length(A) && j <= length(B)) {
49      if (A[i] >= B[j]) {
50         if(A[i-1] >= B[j]) {
51            Tval <- A[i-1]
52            j <- j + 1
53         } else {
54            g <- (B[j] - Tval) / (A[i] - A[i-1])
55            d <- (B[j] - Tval) / (B[j] - B[j-1])
56            pairs <- rbind(pairs, c(g,d))
57            Tval <- B[j]
58            j <- j + 1
59         }
60      } else {
61         if (B[j-1] >= A[i]) {
62            Tval <- B[j-1]
63            i <- i+1
64         } else {
65            g <- (A[i] - Tval) / (A[i] - A[i-1])
66            d <- (A[i] - Tval) / (B[j] - B[j-1])
67            pairs <- rbind(pairs, c(g,d))
68            Tval <- A[i]
69            i <- i + 1
70         }
71      }
72   }
73   # ---------------------
74   # Pairs are computed.
75   # Now compute c_i:
76   # ---------------------
77    mg <- mean(pairs[, 1])
78    md <- mean(pairs[, 2])
79    f1 <- mg*mg+md*md
80    f2 <- sqrt(f1)
81    c_values <- (f1-pairs[, 1]*mg-pairs[, 2]*md)/f2
82    # ----------------------
83    # find  phase function:
84    # ----------------------
85    thetadrift <- mean(abs(c_values))
86    norm_val <- sqrt(2.0)
87    psi <- thetadrift / norm_val
88    #### return phasepairs, c_i, and psi:
89    #list(phasepairs=pairs,
90    #    c_values=c_values,
91    #    psi=psi)
92    #### or just return the phase function value:
93    psi
94 }