请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2726|回复: 4

[经验总结] 最近了解了下Liang-Kleeman信息流

[复制链接]

新浪微博达人勋

发表于 2023-5-20 23:13:13 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
最近了解了下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







密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-5-21 13:17:52 | 显示全部楼层
{:5_213:}{:5_213:}
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2023-5-21 14:23:36 | 显示全部楼层
{:5_213:}{:5_213:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-30 11:50:30 | 显示全部楼层
请问 Liang-Kleeman 输出来的结果 怎么看显著性
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2024-1-3 18:35:34 | 显示全部楼层
本帖最后由 flaggg 于 2024-3-15 16:10 编辑

请问您使用Liang-Kleeman信息流的脚本的时候在pyleo.utils.load_dataset的时候有出现'gbk' codec can't decode byte 0xb0 in position 1703: illegal multibyte sequence
的问题吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表