最近了解了下Liang-Kleeman信息流

1163087557

最近了解了下Liang-Kleeman信息流
发现pyleoclim库中复现了liang教授的代码
import pyleoclim as pyleo

ts_nino=pyleo.utils.load_dataset('NINO3')
ts_air=pyleo.utils.load_dataset('AIR')
pyleo.utils.causality.liang_causality(ts_nino.data, ts_air.data, npt=1, signif_test='isospec', nsim=1000, qs=[0.005, 0.025, 0.05, 0.95, 0.975, 0.995)

即可，同时还自带显著性检验，直接使用还是使用两个一维时间序列。
库里源码
def liang_causality(y1, y2, npt=1, signif_test='isospec', nsim=1000,
                    qs=[0.005, 0.025, 0.05, 0.95, 0.975, 0.995]):
    '''Liang-Kleeman information flow
   
    Estimate the Liang information transfer from series y2 to series y1 with
    significance estimates using either an AR(1) tests with series with the same
    persistence or surrogates with randomized phases.

    Parameters
    ----------

    y1, y2 : array
        vectors of (real) numbers with identical length, no NaNs allowed

    npt : int >=1
        time advance in performing Euler forward differencing,
        e.g., 1, 2. Unless the series are generated with a highly chaotic deterministic system,
        npt=1 should be used
   
    signif_test : str; {'isopersist', 'isospec'}
        the method for significance test
        see signif_isospec and signif_isopersist for details.
        
    nsim : int
        the number of AR(1) surrogates for significance test
        
    qs : list
        the quantiles for significance test

    Returns
    -------

    res : dict
        A dictionary of results including:
            T21 : float
                information flow from y2 to y1 (Note: not y1 -> y2!)
            tau21 : float
                the standardized information flow from y2 to y1
            Z : float
                the total information flow from y2 to y1
            dH1_star : float
                dH*/dt (Liang, 2016)
            dH1_noise : float
            signif_qs :
                the quantiles for significance test
            T21_noise : list
                the quantiles of the information flow from noise2 to noise1 for significance testing
            tau21_noise : list
                the quantiles of the standardized information flow from noise2 to noise1 for significance testing
   
    See also
    --------

    pyleoclim.utils.causality.liang : information flow estimated using the Liang algorithm
   
    pyleoclim.utils.causality.granger_causality : information flow estimated using the Granger algorithm
   
    pyleoclim.utils.causality.signif_isopersist : significance test with AR(1) with same persistence
   
    pyleoclim.utils.causality.causality.signif_isospec : significance test with surrogates with randomized phases
   
    References
    ----------

    Liang, X.S. (2013) The Liang-Kleeman Information Flow: Theory and Applications. Entropy, 15, 327-360, doi:10.3390/e15010327
   
    Liang, X.S. (2014) Unraveling the cause-effect relation between timeseries. Physical review, E 90, 052150
   
    Liang, X.S. (2015) Normalizing the causality between time series. Physical review, E 92, 022126
   
    Liang, X.S. (2016) Information flow and causality as rigorous notions ab initio. Physical review, E 94, 052201

    '''

    dt=1
    nm = np.size(y1)

    grad1 = (y1[0+npt:] - y1[0:-npt]) / (npt)
    grad2 = (y2[0+npt:] - y2[0:-npt]) / (npt)

    y1 = y1[:-npt]
    y2 = y2[:-npt]

    N = nm - npt
    C = np.cov(y1, y2)
    detC = np.linalg.det(C)

    dC = np.ndarray((2, 2))
    dC[0, 0] = np.sum((y1-np.mean(y1))*(grad1-np.mean(grad1)))
    dC[0, 1] = np.sum((y1-np.mean(y1))*(grad2-np.mean(grad2)))
    dC[1, 0] = np.sum((y2-np.mean(y2))*(grad1-np.mean(grad1)))
    dC[1, 1] = np.sum((y2-np.mean(y2))*(grad2-np.mean(grad2)))

    dC /= N-1

    a11 = C[1, 1]*dC[0, 0] - C[0, 1]*dC[1, 0]
    a12 = -C[0, 1]*dC[0, 0] + C[0, 0]*dC[1, 0]

    a11 /= detC
    a12 /= detC

    f1 = np.mean(grad1) - a11*np.mean(y1) - a12*np.mean(y2)
    R1 = grad1 - (f1 + a11*y1 + a12*y2)
    Q1 = np.sum(R1*R1)
    b1 = np.sqrt(Q1*dt/N)

    NI = np.ndarray((4, 4))
    NI[0, 0] = N*dt/b1**2
    NI[1, 1] = dt/b1**2*np.sum(y1*y1)
    NI[2, 2] = dt/b1**2*np.sum(y2*y2)
    NI[3, 3] = 3*dt/b1**4*np.sum(R1*R1) - N/b1**2
    NI[0, 1] = dt/b1**2*np.sum(y1)
    NI[0, 2] = dt/b1**2*np.sum(y2)
    NI[0, 3] = 2*dt/b1**3*np.sum(R1)
    NI[1, 2] = dt/b1**2*np.sum(y1*y2)
    NI[1, 3] = 2*dt/b1**3*np.sum(R1*y1)
    NI[2, 3] = 2*dt/b1**3*np.sum(R1*y2)

    NI[1, 0] = NI[0, 1]
    NI[2, 0] = NI[0, 2]
    NI[2, 1] = NI[1, 2]
    NI[3, 0] = NI[0, 3]
    NI[3, 1] = NI[1, 3]
    NI[3, 2] = NI[2, 3]

    invNI = np.linalg.pinv(NI)
    var_a12 = invNI[2, 2]
    T21 = C[0, 1]/C[0, 0] * (-C[1, 0]*dC[0, 0] + C[0, 0]*dC[1, 0]) / detC
    var_T21 = (C[0, 1]/C[0, 0])**2 * var_a12

    dH1_star= a11
    dH1_noise = b1**2 / (2*C[0, 0])

    Z = np.abs(T21) + np.abs(dH1_star) + np.abs(dH1_noise)

    tau21 = T21 / Z
    dH1_star = dH1_star / Z
    dH1_noise = dH1_noise / Z

    signif_test_func = {
            'isopersist': signif_isopersist,
            'isospec': signif_isospec,
        }

    signif_dict = signif_test_func[signif_test](y1, y2, method='liang', nsim=nsim, qs=qs, npt=npt)
    T21_noise_qs = signif_dict['T21_noise_qs']
    tau21_noise_qs = signif_dict['tau21_noise_qs']

    res = {
        'T21': T21,
        'tau21': tau21,
        'Z': Z,
        'dH1_star': dH1_star,
        'dH1_noise': dH1_noise,
        'signif_qs' : qs,
        'T21_noise' : T21_noise_qs,
        'tau21_noise' : tau21_noise_qs
    }

    return res

def liang(y1, y2, npt=1):
    '''
    Estimate the Liang information transfer from series y2 to series y1

    Parameters
    ----------

    y1, y2 : array
        Vectors of (real) numbers with identical length, no NaNs allowed

    npt : int  >=1
        Time advance in performing Euler forward differencing,
        e.g., 1, 2. Unless the series are generated with a highly chaotic deterministic system,
        npt=1 should be used

    Returns
    -------

    res : dict
        A dictionary of results including:
            T21 : float
                information flow from y2 to y1 (Note: not y1 -> y2!)
            tau21 : float
                the standardized information flow from y2 to y1
            Z : float
                the total information flow from y2 to y1
            dH1_star : float
                dH*/dt (Liang, 2016)
            dH1_noise : float
            
    See also
    --------

    pyleoclim.utils.causality.liang_causality : information flow estimated using the Liang algorithm
    pyleoclim.utils.causality.granger_causality : information flow estimated using the Granger algorithm   
    pyleoclim.utils.causality.signif_isopersist : significance test with AR(1) with same persistence
    pyleoclim.utils.causality.signif_isospec : significance test with surrogates with randomized phases
   
    References
    ----------

    Liang, X.S. (2013) The Liang-Kleeman Information Flow: Theory and
            Applications. Entropy, 15, 327-360, doi:10.3390/e15010327
   
    Liang, X.S. (2014) Unraveling the cause-effect relation between timeseries.
        Physical review, E 90, 052150
   
    Liang, X.S. (2015) Normalizing the causality between time series.
        Physical review, E 92, 022126
   
    Liang, X.S. (2016) Information flow and causality as rigorous notions ab initio.
        Physical review, E 94, 052201

    '''
    dt=1
    nm = np.size(y1)

    grad1 = (y1[0+npt:] - y1[0:-npt]) / (npt)
    grad2 = (y2[0+npt:] - y2[0:-npt]) / (npt)

    y1 = y1[:-npt]
    y2 = y2[:-npt]

    N = nm - npt
    C = np.cov(y1, y2)
    detC = np.linalg.det(C)

    dC = np.ndarray((2, 2))
    dC[0, 0] = np.sum((y1-np.mean(y1))*(grad1-np.mean(grad1)))
    dC[0, 1] = np.sum((y1-np.mean(y1))*(grad2-np.mean(grad2)))
    dC[1, 0] = np.sum((y2-np.mean(y2))*(grad1-np.mean(grad1)))
    dC[1, 1] = np.sum((y2-np.mean(y2))*(grad2-np.mean(grad2)))

    dC /= N-1

    a11 = C[1, 1]*dC[0, 0] - C[0, 1]*dC[1, 0]
    a12 = -C[0, 1]*dC[0, 0] + C[0, 0]*dC[1, 0]

    a11 /= detC
    a12 /= detC

    f1 = np.mean(grad1) - a11*np.mean(y1) - a12*np.mean(y2)
    R1 = grad1 - (f1 + a11*y1 + a12*y2)
    Q1 = np.sum(R1*R1)
    b1 = np.sqrt(Q1*dt/N)

    NI = np.ndarray((4, 4))
    NI[0, 0] = N*dt/b1**2
    NI[1, 1] = dt/b1**2*np.sum(y1*y1)
    NI[2, 2] = dt/b1**2*np.sum(y2*y2)
    NI[3, 3] = 3*dt/b1**4*np.sum(R1*R1) - N/b1**2
    NI[0, 1] = dt/b1**2*np.sum(y1)
    NI[0, 2] = dt/b1**2*np.sum(y2)
    NI[0, 3] = 2*dt/b1**3*np.sum(R1)
    NI[1, 2] = dt/b1**2*np.sum(y1*y2)
    NI[1, 3] = 2*dt/b1**3*np.sum(R1*y1)
    NI[2, 3] = 2*dt/b1**3*np.sum(R1*y2)

    NI[1, 0] = NI[0, 1]
    NI[2, 0] = NI[0, 2]
    NI[2, 1] = NI[1, 2]
    NI[3, 0] = NI[0, 3]
    NI[3, 1] = NI[1, 3]
    NI[3, 2] = NI[2, 3]

    invNI = np.linalg.pinv(NI)
    var_a12 = invNI[2, 2]
    T21 = C[0, 1]/C[0, 0] * (-C[1, 0]*dC[0, 0] + C[0, 0]*dC[1, 0]) / detC
    var_T21 = (C[0, 1]/C[0, 0])**2 * var_a12

    dH1_star= a11
    dH1_noise = b1**2 / (2*C[0, 0])

    Z = np.abs(T21) + np.abs(dH1_star) + np.abs(dH1_noise)

    tau21 = T21 / Z
    dH1_star = dH1_star / Z
    dH1_noise = dH1_noise / Z

    res = {
        'T21': T21,
        'tau21': tau21,
        'Z': Z,
        'dH1_star': dH1_star,
        'dH1_noise': dH1_noise,
    }

    return res

一大碗年糕

{:5_213:}{:5_213:}

内马尔

{:5_213:}{:5_213:}

王启璐

请问 Liang-Kleeman 输出来的结果 怎么看显著性

flaggg

请问您使用Liang-Kleeman信息流的脚本的时候在pyleo.utils.load_dataset的时候有出现'gbk' codec can't decode byte 0xb0 in position 1703: illegal multibyte sequence
的问题吗？