Return to Question

UPDATE: The solution I tried from the answers below improved runtime for larger inputs than my original approach. However, I still face the problem of running this statistical function a 1000 times per array pair. The background for this is to perform a permutation test on the statistic with 'x' (1000) number of permutations to arrive at a significance level at p= .01 Some asked below why I need to run this function billions of times. So my data is a timeseries data for events: 'E1,E2,E3.....En' I need to find combinations of all possible event pairs from this data i.e. nCr. (E1,E2), (E2,E3)....nCr times. This leads to a large number of event pairs about billions of event pairs for my data. Add to that the time complexity of a permutation test (1000x times/pair). Any suggestions if any other parts like permutation tests can be improved. Thanks

UPDATE: The solution I tried from the answers below improved runtime for larger inputs than my original approach. However, I still face the problem of running this statistical function a 1000 times per array pair. The background for this is to perform a permutation test on the statistic with 'x' (1000) number of permutations to arrive at a significance level at p= .01 Some asked below why I need to run this function billions of times. So my data is a timeseries data for events: 'E1,E2,E3.....En' I need to find combinations of all possible event pairs from this data i.e. nCr. (E1,E2), (E2,E3)....nCr times. This leads to a large number of event pairs about billions of event pairs for my data. Add to that the time complexity of a permutation test (1000x times/pair). Any suggestions if any other parts like permutation tests can be improved. Thanks
UPDATE2: Is size always ~8000? No it depends on input data, array size can be between 2000 ~ 10,000. What's your tau value? tau is a random positive integer or float and we jus test for different tau values, so I cannot given a certain Tau.

#e.g. inputs
#ts_1 = [0,1,1,0,0,1,1,0,......] #timeseries1
#ts_2 = [1,1,1,1,1,1,1,0,......] #timeseries2
# tau = positive integer/float   #time lag

import numpy as np
import pandas as pd
from itertools import product

 def eps(ts_1, ts_2, tau):  
    Q_tau = 0
    q_tau = 0
    
    event_index1, = np.where(np.array(ts_1) == 0)
    n1 = event_index1.shape[0]
    
    event_index2, = np.where(np.array(ts_2) == 0)
    n2 = event_index2.shape[0]
    
    if (n1 != 0 and n2 != 0):
        matching_idx = set(event_index1).intersection(event_index2)
        c_ij = c_ji = 0.5 *len(matching_idx)
        
        for x,y in product(event_index1,event_index2):
            if x-y > 0 and (x-y)<= tau:
                c_ij += 1
            elif y-x > 0 and (y-x) <= tau:
                c_ji += 1                     
                  
        Q_tau = (c_ij+c_ji)/math.sqrt( n1 * n2 )
        q_tau = (c_ij - c_ji)/math.sqrt( n1 * n2 )
    
    return Q_tau, q_tau

#permutation test

 def eps_permutationtest(ts_1,ts_2,permutations=1000,Q_tau): #Q_tau= original statistic to evaluate significance for
        
        lst_qij = []    # list to hold Q_tau values from each shuffle 
        
        for i in range(0,permutations):
            np.random.shuffle(ts_1)         #shuffle ts_1 everytime
            q_ij = eps(ts_1 = ts_1, ts_2 = ts_2 , tau= 5) 
            lst_qij.append(q_ij)
        lst_qij = np.array(lst_qij)
        
    
        p_val = len(np.where(lst_qij >= [Q_tau])[0])/permutations
        
        
        return p_val  # significant if p < .01 alpha

#e.g. inputs
#ts_1 = [0,1,1,0,0,1,1,0,......] #timeseries1
#ts_2 = [1,1,1,1,1,1,1,0,......] #timeseries2
# tau = positive integer/float   #time lag

import numpy as np
from itertools import product

 def eps(ts_1, ts_2, tau):  
    Q_tau = 0
    q_tau = 0
    
    event_index1, = np.where(np.array(ts_1) == 0)
    n1 = event_index1.shape[0]
    
    event_index2, = np.where(np.array(ts_2) == 0)
    n2 = event_index2.shape[0]
    
    if (n1 != 0 and n2 != 0):
        matching_idx = set(event_index1).intersection(event_index2)
        c_ij = c_ji = 0.5 *len(matching_idx)
        
        for x,y in product(event_index1,event_index2):
            if x-y > 0 and (x-y)<= tau:
                c_ij += 1
            elif y-x > 0 and (y-x) <= tau:
                c_ji += 1                     
                  
        Q_tau = (c_ij+c_ji)/math.sqrt( n1 * n2 )
        q_tau = (c_ij - c_ji)/math.sqrt( n1 * n2 )
    
    return Q_tau, q_tau

#permutation test

 def eps_permutationtest(ts_1,ts_2,permutations=1000,Q_tau): #Q_tau= original statistic to evaluate significance for
        
     lst_Q_tau = []    # list to hold Q_tau values from each shuffle 
        
     for i in range(0,permutations):
         np.random.shuffle(ts_1)         #shuffle ts_1 everytime
         Q_tau = eps(ts_1 = ts_1, ts_2 = ts_2 , tau= 5) 
         lst_Q_tau.append(Q_tau)
     lst_Q_tau= np.array(lst_Q_tau)
        
     p_val = len(np.where(lst_Q_tau>= [Q_tau])[0])/permutations
        
        
     return p_val  # significant if p < .01 alpha

I write a function which takes as input 2 arrays of zeros and ones ~8000 elements per array. Now the array size is not fixed it depends on the input data and can be much higher number of elements. The density of ones is expected to be higher than zeros for most cases, but outliers may exist. My function eps calculates a statistic on these arrays and returns the output. This statistic is just a way to find if 2 timeseries correlate with each other. Relatively trivial operations of just checking for 0 and noting the index where 0 is found in array and then some calculations. I tried my best to optimize for speed but the best I could get is 4.5 ~5 seconds (for 18k array pairs) using timeit (1000 runs) library. Fast execution is important as I need to run this function on billions of array pairs.

UPDATE: The solution I tried from the answers below improved well for larger inputs than the original approach. However, I still face the problem of running this statistical function a 1000 times per array pair. The background to this is to perform a permutation test on the statistic with 'x' (1000) number of permutations to arrive get a significance level at p= .01 Some asked below why I need to run this function billions of times. So my data is a timeseries data for an event 'E1,E2,E3.....En' I need to find combinations of 2 event pairs from this data i.e. nCr. 'E1,E2', 'E2,E3'....nCr times. This leading to a large number of event pairs about billions of event pairs for my data. Add to that the time complexity of a permutation test (1000x times). Any suggestions if any other parts like permutation test can be improved. Thanks

#e.g. inputs
#ts_1 = [0,1,1,0,0,1,1,0,......]
#ts_2 = [1,1,1,1,1,1,1,0,......]
# tau = positive integer/float

import numpy as np
import pandas as pd
from itertools import product

 def eps(ts_1, ts_2, tau):  
    Q_tau = 0
    q_tau = 0
    
    event_index1, = np.where(np.array(ts_1) == 0)
    n1 = event_index1.shape[0]
    
    event_index2, = np.where(np.array(ts_2) == 0)
    n2 = event_index2.shape[0]
    
    if (n1 != 0 and n2 != 0):
        matching_idx = set(event_index1).intersection(event_index2)
        c_ij = c_ji = 0.5 *len(matching_idx)
        
        for x,y in product(event_index1,event_index2):
            if x-y > 0 and (x-y)<= tau:
                c_ij += 1
            elif y-x > 0 and (y-x) <= tau:
                c_ji += 1                     
                  
        Q_tau = (c_ij+c_ji)/math.sqrt( n1 * n2 )
        q_tau = (c_ij - c_ji)/math.sqrt( n1 * n2 )
    
    return Q_tau, q_tau

I write a function which takes as input 2 arrays of zeros and ones ~8000 elements per array. Input array size is variable, can have much higher number of elements. The density of ones is expected to be higher than zeros for most cases, but outliers may exist. My function eps calculates a statistic on these array pairs and returns the output. This statistic is just a way to find if 2 timeseries data correlate with each other. Relatively trivial operations of checking for 0 and noting the index where 0 is found in array and then some calculations. I tried my best to optimize for speed but the best I could get is 4.5 ~5 seconds (for 18k array pairs) using timeit (1000 runs) library. Fast execution is important as I need to run this function on billions of array pairs.

UPDATE: The solution I tried from the answers below improved runtime for larger inputs than my original approach. However, I still face the problem of running this statistical function a 1000 times per array pair. The background for this is to perform a permutation test on the statistic with 'x' (1000) number of permutations to arrive at a significance level at p= .01 Some asked below why I need to run this function billions of times. So my data is a timeseries data for events: 'E1,E2,E3.....En' I need to find combinations of all possible event pairs from this data i.e. nCr. (E1,E2), (E2,E3)....nCr times. This leads to a large number of event pairs about billions of event pairs for my data. Add to that the time complexity of a permutation test (1000x times/pair). Any suggestions if any other parts like permutation tests can be improved. Thanks

#e.g. inputs
#ts_1 = [0,1,1,0,0,1,1,0,......] #timeseries1
#ts_2 = [1,1,1,1,1,1,1,0,......] #timeseries2
# tau = positive integer/float   #time lag

import numpy as np
import pandas as pd
from itertools import product

 def eps(ts_1, ts_2, tau):  
    Q_tau = 0
    q_tau = 0
    
    event_index1, = np.where(np.array(ts_1) == 0)
    n1 = event_index1.shape[0]
    
    event_index2, = np.where(np.array(ts_2) == 0)
    n2 = event_index2.shape[0]
    
    if (n1 != 0 and n2 != 0):
        matching_idx = set(event_index1).intersection(event_index2)
        c_ij = c_ji = 0.5 *len(matching_idx)
        
        for x,y in product(event_index1,event_index2):
            if x-y > 0 and (x-y)<= tau:
                c_ij += 1
            elif y-x > 0 and (y-x) <= tau:
                c_ji += 1                     
                  
        Q_tau = (c_ij+c_ji)/math.sqrt( n1 * n2 )
        q_tau = (c_ij - c_ji)/math.sqrt( n1 * n2 )
    
    return Q_tau, q_tau

#permutation test

 def eps_permutationtest(ts_1,ts_2,permutations=1000,Q_tau): #Q_tau= original statistic to evaluate significance for
        
        lst_qij = []    # list to hold Q_tau values from each shuffle 
        
        for i in range(0,permutations):
            np.random.shuffle(ts_1)         #shuffle ts_1 everytime
            q_ij = eps(ts_1 = ts_1, ts_2 = ts_2 , tau= 5) 
            lst_qij.append(q_ij)
        lst_qij = np.array(lst_qij)
        
    
        p_val = len(np.where(lst_qij >= [Q_tau])[0])/permutations
        
        
        return p_val  # significant if p < .01 alpha