<jupyter_start><jupyter_text>Maximum Likelihood Estimation [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/danhtaihoang/network-inference/master?filepath=sphinx%2Fcodesource%2Fmle.ipynb)Maximum likelihood estimation (MLE) is a standard method in the field of network inference. The algorithm of MLE for the kinetic model contains the following steps:(i) Initialize $W_{ij}$ at random;(ii) Compute local field $H_i(t) = \sum_j W_{ij} \sigma_j (t)$;(iii) Compute probability;$$P[\sigma_i(t+1)|\vec{\sigma}(t)] = \frac{\exp [ \sigma_i(t+1) H_i(t)]}{\sum_{\sigma_i(t+1)} \exp[\sigma_i(t+1) H_i(t)]}$$(iv) Compute likelihood$\mathcal{P}_i = \prod_{t=1}^{L-1}P[\sigma_i(t+1)|\vec{\sigma}(t)]$;(v) Update $W_{ij}$ incrementally through$W_{ij}^{\textrm{new}} = W_{ij} + \frac{\alpha}{L-1} \frac{\partial \log {\cal{P}}_i }{\partial W_{ij}}$where $\alpha$ represent learning rate, and$\frac{\partial \log {\mathcal{P}}_i}{\partial W_{ij}} = \sum_{t=1}^{L-1} \lbrace \sigma_i(t+1) \sigma_j(t) - \tanh [H_i(t)] \sigma_j(t) \rbrace$;(vi) Repeat (ii)-(v) for some iterations;(vii) Perform (ii)-(vi) in parallel for every variable index$i \in \{1, 2, \cdots, N\}$.In the following, we will test the performance of this method for various cases, small coupling variability and large coupling variability. We will show that the MLE becomes over fitting in the regime of large coupling variability and small sample sizes. Applying our stopping criterion to MLE can make this method works significantly faster and slighlty more accurate.First of all, we import the packages to the notebook:<jupyter_code>import numpy as np
import sys
import timeit
import matplotlib.pyplot as plt
import simulate
import inference
%matplotlib inline

np.random.seed(1)<jupyter_output><empty_output><jupyter_text>We first consider a system of $N = 100$ variables. Data length $L=2000$ is used for the illustration. Leaning rate $\alpha=1$ is selected.<jupyter_code># parameter setting:
n = 100
l = 2000
rate = 1.<jupyter_output><empty_output><jupyter_text>Small coupling variabilityFor the regime of small coupling variability, as other methods, we select $g = 2.0$ for the demonstration.<jupyter_code>g = 2.0<jupyter_output><empty_output><jupyter_text>The coupling matrix is generated acorrding to the parameters.<jupyter_code>w0 = np.random.normal(0.0,g/np.sqrt(n),size=(n,n))<jupyter_output><empty_output><jupyter_text>From the coupling matrix, we generate time series of variables based on the kinetic ising model:<jupyter_code>s = simulate.generate_data(w0,l)<jupyter_output><empty_output><jupyter_text>Similar to FEM and eMF, MLE is an iterative algorithm. We first run this method as discribed in the literature, without stopping criterion.<jupyter_code>start_time = timeit.default_timer()

w = inference.mle(s,rate,stop_criterion='no')

stop_time=timeit.default_timer()
run_time=stop_time-start_time
print('run_time:',run_time)<jupyter_output>('run_time:', 158.09815406799316)<jupyter_text>The inferred couplings are plotted to compare with actual couplings:<jupyter_code>plt.figure(figsize=(11,3.2))

plt.subplot2grid((1,3),(0,0))
plt.title('actual coupling matrix')
plt.imshow(w0,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,1))
plt.title('predicted coupling matrix')
plt.imshow(w,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,2))
plt.title('predicted couplings vs. actual couplings')
plt.plot([-1,1],[-1,1],'r--')
plt.scatter(w0,w)
plt.xlabel('actual couplings')
plt.ylabel('inferred couplings')

plt.tight_layout(h_pad=1, w_pad=1.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Mean square error is calculated as the following:<jupyter_code>MSE = np.mean((w-w0)**2)
print('MSE:',MSE)<jupyter_output>('MSE:', 0.0030760425560816744)<jupyter_text>Now, we apply our stopping criterion for this method. According to this criterion, we stop the iteration when the discrepancy between observed variables and its model expectations starts to increase.<jupyter_code>start_time = timeit.default_timer()

w = inference.mle(s,rate,stop_criterion='yes')

stop_time=timeit.default_timer()
run_time=stop_time-start_time
print('run_time:',run_time)<jupyter_output>('run_time:', 8.734601020812988)<jupyter_text>The inferred couplings are plotted:<jupyter_code>plt.figure(figsize=(11,3.2))

plt.subplot2grid((1,3),(0,0))
plt.title('actual coupling matrix')
plt.imshow(w0,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,1))
plt.title('predicted coupling matrix')
plt.imshow(w,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,2))
plt.title('predicted couplings vs. actual couplings')
plt.plot([-1,1],[-1,1],'r--')
plt.scatter(w0,w)
plt.xlabel('actual couplings')
plt.ylabel('inferred couplings')

plt.tight_layout(h_pad=1, w_pad=1.5)
plt.show()
MSE = np.mean((w-w0)**2)
print('MSE:',MSE)<jupyter_output>('MSE:', 0.0029889738951087916)<jupyter_text>So, our stopping criterion makes MLE works slightly better (in accuracy) and significantly faster.  Large coupling variabilityNow, let us consider the regime of large coupling variability, $g=4$ for instance.<jupyter_code>g = 4.0
w0 = np.random.normal(0.0,g/np.sqrt(n),size=(n,n))

s = simulate.generate_data(w0,l)<jupyter_output><empty_output><jupyter_text>When our stopping criterion does not apply to MLE:<jupyter_code>start_time = timeit.default_timer()

w = inference.mle(s,rate,stop_criterion='no')

stop_time=timeit.default_timer()
run_time=stop_time-start_time
print('run_time:',run_time)
plt.figure(figsize=(11,3.2))

plt.subplot2grid((1,3),(0,0))
plt.title('actual coupling matrix')
plt.imshow(w0,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,1))
plt.title('predicted coupling matrix')
plt.imshow(w,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,2))
plt.title('predicted couplings vs. actual couplings')
plt.plot([-2.0,2.0],[-2.0,2.0],'r--')
plt.scatter(w0,w)
plt.xlabel('actual couplings')
plt.ylabel('inferred couplings')

plt.tight_layout(h_pad=1, w_pad=1.5)
plt.show()
MSE = np.mean((w-w0)**2)
print('MSE:',MSE)<jupyter_output>('MSE:', 0.01781885966000276)<jupyter_text>When our stopping criterion does apply to MLE:<jupyter_code>start_time = timeit.default_timer()

w = inference.mle(s,rate,stop_criterion='yes')

stop_time=timeit.default_timer()
run_time=stop_time-start_time
print('run_time:',run_time)
plt.figure(figsize=(11,3.2))

plt.subplot2grid((1,3),(0,0))
plt.title('actual coupling matrix')
plt.imshow(w0,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,1))
plt.title('predicted coupling matrix')
plt.imshow(w,cmap='rainbow',origin='lower')
plt.xlabel('j')
plt.ylabel('i')
plt.clim(-0.5,0.5)
plt.colorbar(fraction=0.045, pad=0.05,ticks=[-0.5,0,0.5])

plt.subplot2grid((1,3),(0,2))
plt.title('predicted couplings vs. actual couplings')
plt.plot([-2.0,2.0],[-2.0,2.0],'r--')
plt.scatter(w0,w)
plt.xlabel('actual couplings')
plt.ylabel('inferred couplings')

plt.tight_layout(h_pad=1, w_pad=1.5)
plt.show()
MSE = np.mean((w-w0)**2)
print('MSE:',MSE)<jupyter_output>('MSE:', 0.0172333170396427)
{'fraction_non_alphanumeric': 0.18089613557581205, 'fraction_numerical': 0.05289510848632687, 'mean_word_length': 7.340471092077088, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '250349', 'license': 'MIT', 'n_tokens_mistral': 3400, 'n_tokens_neox': 3048, 'n_words': 555, 'path': 'sphinx/codesource/.ipynb_checkpoints/mle-checkpoint.ipynb', 'repo_name': 'danhtaihoang/network-inference', 'size': '7789'}