R function
- phasefunc.R
To make the function available for R environment:
> source('phasefunc.R')
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 }