causality_for_spaceweather/notebooks/_build/jupyter_execute/chapter1_it/tigramite_basic.py

#!/usr/bin/env python
# coding: utf-8

# # Tigramite基本使用教程

# ## 线性测试

# In[1]:


import numpy as np
import matplotlib
from matplotlib import pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
## use `%matplotlib notebook` for interactive figures
# plt.style.use('ggplot')
import sklearn

import tigramite
from tigramite import data_processing as pp
from tigramite import plotting as tp
from tigramite.pcmci import PCMCI
from tigramite.independence_tests import ParCorr, GPDC, CMIknn, CMIsymb


# 数据生成函数：
# 
# \begin{aligned}
# X^0_t &= 0.7 X^0_{t-1} - 0.8 X^1_{t-1} + \eta^0_t\\
# X^1_t &= 0.8 X^1_{t-1} + 0.8 X^3_{t-1} + \eta^1_t\\
# X^2_t &= 0.5 X^2_{t-1} + 0.5 X^1_{t-2} + 0.6 X^3_{t-3} + \eta^2_t\\
# X^3_t &= 0.7 X^3_{t-1} + \eta^3_t\\
# \end{aligned}
# 
# 其中$\eta$满足均值为0方差为1的高斯分布，是添加的随机噪音。
# 
# 目标：
# 
# 重构出每一个变量的`drivers`。
# 

# In[2]:


np.random.seed(42)     
links_coeffs = {0: [((0, -1), 0.7), ((1, -1), -0.8)],
                1: [((1, -1), 0.8), ((3, -1), 0.8)],
                2: [((2, -1), 0.5), ((1, -2), 0.5), ((3, -3), 0.6)],
                3: [((3, -1), 0.4)],
                }
T = 1000    
data, true_parents_neighbors = pp.var_process(links_coeffs, T=T)


# In[3]:


T, N = data.shape
var_names = [r'$X^0$', r'$X^1$', r'$X^2$', r'$X^3$']
dataframe = pp.DataFrame(data, 
                         datatime = np.arange(len(data)), 
                         var_names=var_names)


# In[4]:


data.shape


# In[5]:


tp.plot_timeseries(dataframe)
plt.show()


# In[6]:


parcorr = ParCorr(significance='analytic')
pcmci = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=parcorr,
    verbosity=1)


# In[7]:


correlations = pcmci.get_lagged_dependencies(tau_max=20, val_only=True)['val_matrix']
lag_func_matrix = tp.plot_lagfuncs(val_matrix=correlations,
                                   setup_args={'figsize': (6,6),'var_names':var_names,
                                    'x_base':10, 'y_base':.5})


# In[8]:


pcmci.verbosity = 0
results = pcmci.run_pcmci(tau_max=8, pc_alpha=None)


# In[9]:


print("p-values")
print (results['p_matrix'].round(3))
print("MCI partial correlations")
print (results['val_matrix'].round(2))


# In[10]:


q_matrix = pcmci.get_corrected_pvalues(p_matrix=results['p_matrix'], fdr_method='fdr_bh')
pcmci.print_significant_links(
        p_matrix = results['p_matrix'], 
        q_matrix = q_matrix,
        val_matrix = results['val_matrix'],
        alpha_level = 0.01)


# In[11]:


link_matrix = pcmci.return_significant_links(pq_matrix=q_matrix,
                        val_matrix=results['val_matrix'], alpha_level=0.01)['link_matrix']


# In[12]:


link_matrix.shape


# In[13]:


tp.plot_graph(
    figsize=(5,5),
    val_matrix=results['val_matrix'],
    link_matrix=link_matrix,
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    ); plt.show()


# In[14]:


# Plot time series graph
tp.plot_time_series_graph(
    val_matrix=results['val_matrix'],
    link_matrix=link_matrix,
    var_names=var_names,
    link_colorbar_label='MCI',
    ); plt.show()


# ## 非线性测试

# In[15]:


import numpy as np
import matplotlib
from matplotlib import pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
## use `%matplotlib notebook` for interactive figures
# plt.style.use('ggplot')
import sklearn

import tigramite
from tigramite import data_processing as pp
from tigramite import plotting as tp
from tigramite.pcmci import PCMCI
from tigramite.independence_tests import ParCorr, GPDC, CMIknn, CMIsymb


# 
# \begin{align*}
#     X^0_t &= 0.2 (X^1_{t-1})^2 + \eta^0_t\\
#     X^1_t &= \eta^1_t \\
#     X^2_t &= 0.3 (X^1_{t-2})^2 + \eta^2_t
# \end{align*}

# In[16]:


var_names = [r'$X^0$', r'$X^1$', r'$X^2$']


# In[17]:


np.random.seed(1)
data = np.random.randn(500, 3)
for t in range(1, 500):
    data[t, 0] += 0.4*data[t-1, 1]**2
    data[t, 2] += 0.3*data[t-2, 1]**2
dataframe = pp.DataFrame(data, var_names=var_names)
tp.plot_timeseries(dataframe)
plt.show()


# ### corr

# In[18]:


parcorr = ParCorr(significance='analytic')
pcmci_parcorr = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=parcorr,
    verbosity=0)


# In[19]:


correlations = pcmci_parcorr.get_lagged_dependencies(tau_max=20, val_only=True)['val_matrix']


# In[20]:


lag_func_matrix = tp.plot_lagfuncs(val_matrix=correlations,
                                   setup_args={'var_names':var_names, 'figsize': (5,5),
                                    'x_base':5, 'y_base':.5})


# In[21]:


results = pcmci_parcorr.run_pcmci(tau_max=2, pc_alpha=0.2)
pcmci_parcorr.print_significant_links(
        p_matrix = results['p_matrix'], 
        val_matrix = results['val_matrix'],
        alpha_level = 0.01)


# In[22]:


q_matrix = pcmci_parcorr.get_corrected_pvalues(p_matrix=results['p_matrix'], fdr_method='fdr_bh')

link_matrix = pcmci_parcorr.return_significant_links(pq_matrix=q_matrix,
                        val_matrix=results['val_matrix'], alpha_level=0.01)['link_matrix']


# In[23]:


tp.plot_graph(
    val_matrix=results['val_matrix'],
    link_matrix=link_matrix,
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    ); plt.show()


# ### gpdc

# In[24]:


gpdc = GPDC(significance='analytic', gp_params=None)
pcmci_gpdc = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=gpdc,
    verbosity=0)


# In[25]:


results = pcmci_gpdc.run_pcmci(tau_max=2, pc_alpha=0.1)
pcmci_gpdc.print_significant_links(
        p_matrix = results['p_matrix'], 
        val_matrix = results['val_matrix'],
        alpha_level = 0.01)


# In[26]:


q_matrix = pcmci_parcorr.get_corrected_pvalues(p_matrix=results['p_matrix'], fdr_method='fdr_bh')

link_matrix = pcmci_parcorr.return_significant_links(pq_matrix=q_matrix,
                        val_matrix=results['val_matrix'], alpha_level=0.01)['link_matrix']


# In[27]:


tp.plot_graph(
    val_matrix=results['val_matrix'],
    link_matrix=link_matrix,
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    ); plt.show()


# ### cmiknn

# In[28]:


cmi_knn = CMIknn(significance='shuffle_test', knn=0.1, shuffle_neighbors=5, transform='ranks')
pcmci_cmi_knn = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=cmi_knn,
    verbosity=2)
results = pcmci_cmi_knn.run_pcmci(tau_max=2, pc_alpha=0.05)
pcmci_cmi_knn.print_significant_links(
        p_matrix = results['p_matrix'], 
        val_matrix = results['val_matrix'],
        alpha_level = 0.01)


# In[29]:


q_matrix = pcmci_parcorr.get_corrected_pvalues(p_matrix=results['p_matrix'], fdr_method='fdr_bh')

link_matrix = pcmci_parcorr.return_significant_links(pq_matrix=q_matrix,
                        val_matrix=results['val_matrix'], alpha_level=0.01)['link_matrix']


# In[30]:


tp.plot_graph(
    val_matrix=results['val_matrix'],
    link_matrix=link_matrix,
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    ); plt.show()


# In[ ]:





# In[ ]: