In [34]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from IPython.display import Image

%pylab inline

pylab.rcParams['figure.figsize'] = (10, 6)

%load_ext autoreload
%autoreload 2
Populating the interactive namespace from numpy and matplotlib
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

1. Install TISEAN

2. The standard first step in the analysis of a scalar time series from a nonlinear dynamical system is to choose the delay $\tau$. Use TISEAN's mutual tool to construct a plot of mutual information vesus $\tau$ for data2.first250sec. We're looking for the first minumum.

First we run TISEAN with default parameters.

In [35]:
%%bash --out p2_mutual_string --err error

mutual ps8data/data2.first250sec

Using iPython bash magic we export the output to stdout (there are only 20 rows) and import into a pandas DataFrame.

In [36]:
p2_mutual = np.array([r.split(' ') for r in
                p2_mutual_string.split('\n')][1:-1],
                dtype=np.float64)

p2d = pd.DataFrame(p2_mutual, columns=['Tau', 'AMI'])

p2d.plot(x='Tau', y='AMI', label=r'AMI vs. $\tau$', legend=False)
plt.show()

We don't see a minimum, so we're going to slowly increase the max time delay parameter until we see that minimum.

In [37]:
%%bash --out p2_mutual_string --err error

mutual -D 200 ps8data/data2.first250sec
In [38]:
p2_mutual = np.array([r.split(' ') for r in
                p2_mutual_string.split('\n')][1:-1],
                dtype=np.float64)

p2d = pd.DataFrame(p2_mutual, columns=['Tau', 'AMI'])

p2d.plot(x='Tau', y='AMI', label=r'AMI vs. $\tau$', legend=False)
plt.show()

This minimum is now apparent, so let's figure out the $\tau$ value at which it occurs.

In [39]:
p2min = p2d.idxmin()
p2min
Out[39]:
Tau      0
AMI    155
dtype: int64

We need the data time interval to figure out this value in seconds.

In [40]:
p2_file_data = pd.read_csv('ps8data/data2.first250sec',
                            names=['angle', '', '', 'time'],
                            sep=' ')[['angle', 'time']]
interval = p2_file_data['time'][1] - p2_file_data['time'][0]
interval
Out[40]:
0.002

Yielding a $\tau$ value of

In [41]:
p2min['AMI'] * interval
Out[41]:
0.31

3. The standard second step in nonlinear time-series analysis is to choose the embedding dimension $m$. Use TISEAN's false_nearest tool to construct a plot of the percent of false-near neighbors versus $m$ for the same dataset. Use an $m$ range of $[1, 10]$, using $\tau$ from before, and leave the other parameters at their default values. Plot the results.

The first step is to use TISEAN.

In [42]:
%%bash --out p3_false_nearest_string --err error

false_nearest -d 155 -m 1 -M 1,10 ps8data/data2.first250sec
In [43]:
p3d = pd.DataFrame(np.array([r.split(' ') for r in
                   p3_false_nearest_string.split('\n')[:-1]],
                   dtype=np.float64),
                   columns=['Dimension', 'Frac. FNN',
                          'Avg. Neighborhood Size',
                          'Avg. Neighborhood^2 Size'])
p3d
Out[43]:
Dimension Frac. FNN Avg. Neighborhood Size Avg. Neighborhood^2 Size
0 1 0.893131 0.001500 0.001500
1 2 0.659106 0.001569 0.003657
2 3 0.475366 0.001912 0.004873
3 4 0.335245 0.002482 0.005668
4 5 0.232490 0.003070 0.006467
5 6 0.163109 0.003606 0.007233
6 7 0.119633 0.004089 0.007943
7 8 0.090941 0.004534 0.008621
8 9 0.070794 0.004944 0.009264
9 10 0.056035 0.005338 0.009888

Now we can plot this.

In [44]:
p3d.plot(x='Dimension', y='Frac. FNN',
         label='False Nearest Neighbors', legend=False)
plt.plot(np.arange(0, 11), 0.1 * np.ones(11), 'k--')
plt.show()

Looking at the data we see that $m=8$ is the optimal value.

4. Use TISEAN's delay tool to embed the same dataset using the found values of $\tau$ and $m$. Compare to the plot from PS8.

First with TISEAN.

In [45]:
%%bash --out p4_delay_string --err error

delay -m 8 -d 155 ps8data/data2.first250sec
In [46]:
p4d = pd.DataFrame(np.array([r.split(' ')[:-1] for r in
                   p4_delay_string.split('\n')[:-1]],
                   dtype=np.float64))
plt.figure(figsize=(50, 50))
plt.scatter(p4d[0], p4d[1])
plt.show()