{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Introduction\n### This is a follow-up kernel of [my previous LANL EDA kernel](https://www.kaggle.com/tarunpaparaju/lanl-earthquake-prediction-fresh-eda). In this kernel, I explore even more new features and visualize their relationships with the target.\n"},{"metadata":{},"cell_type":"markdown","source":"### Please upvote this kernel if you like it :)"},{"metadata":{},"cell_type":"markdown","source":"<center><img src=\"https://i.imgur.com/hBPv3fh.png\" width=\"750px\"></center>"},{"metadata":{},"cell_type":"markdown","source":"### Import necessary libraries"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import os\nimport gc\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\ntqdm.pandas()\n\nfrom numba import jit\nfrom math import log, floor\nfrom sklearn.neighbors import KDTree\nfrom scipy.signal import periodogram, welch\n\nfrom keras.layers import *\nfrom keras.models import *\nfrom tqdm import tqdm\nfrom sklearn.model_selection import train_test_split \nfrom keras import backend as K\nfrom keras import optimizers\nfrom sklearn.model_selection import GridSearchCV, KFold\nfrom keras.callbacks import *\nfrom keras import activations\nfrom keras import regularizers\nfrom keras import initializers\nfrom keras import constraints\nfrom keras.engine import Layer\nfrom keras.engine import InputSpec\nfrom keras.objectives import categorical_crossentropy\nfrom keras.objectives import sparse_categorical_crossentropy\nfrom keras.utils import plot_model\nfrom keras.utils.vis_utils import model_to_dot\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom IPython.display import SVG\n\nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Initialize necessay constants"},{"metadata":{"trusted":true},"cell_type":"code","source":"SIGNAL_LEN = 150000\nMIN_NUM = -27\nMAX_NUM = 28","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Download seismic signal data along with targets (time left for occurance of laboratory earthquake)"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"seismic_signals = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Extract the acoustic data and targets from the dataframe\nNote : I delete the original dataframe to save memory"},{"metadata":{"trusted":true},"cell_type":"code","source":"acoustic_data = seismic_signals.acoustic_data\ntime_to_failure = seismic_signals.time_to_failure\ndata_len = len(seismic_signals)\ndel seismic_signals\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Break the data down into parts\nWe have one long array of seismic data. We will break it down into chunks of size 150k (**SIGNAL_LEN**) and each chunk will be one signal in our data (this is because each segment in the test data has length 150k). The **time_to_failure** at the last time step of each segment becomes the target associated with that segment."},{"metadata":{"trusted":true},"cell_type":"code","source":"signals = []\ntargets = []\n\nfor i in range(data_len//SIGNAL_LEN):\n    min_lim = SIGNAL_LEN * i\n    max_lim = min([SIGNAL_LEN * (i + 1), data_len])\n    \n    signals.append(list(acoustic_data[min_lim : max_lim]))\n    targets.append(time_to_failure[max_lim])\n    \ndel acoustic_data\ndel time_to_failure\ngc.collect()\n    \nsignals = np.array(signals)\ntargets = np.array(targets)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Functions for preparing signal features"},{"metadata":{},"cell_type":"markdown","source":"### Scaling the signals\nThis function scales the seismic signals from its original range (-27 to 28 : this where 99% of the data lies) to the range (-1, 1)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def min_max_transfer(ts, min_value, max_value, range_needed=(-1,1)):\n    ts_std = (ts - min_value) / (max_value - min_value)\n\n    if range_needed[0] < 0:    \n        return ts_std * (range_needed[1] + abs(range_needed[0])) + range_needed[0]\n    else:\n        return ts_std * (range_needed[1] - range_needed[0]) + range_needed[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Extracting features from each part of the segment\nThe original long seismic signal has already been broken down into several segments. The segments are scaled using the **min_max_transfer** function. Then, we break down each segment into several parts. Usual features such as mean, standard deviation, range, percentiles etc are calculated over each part of the segment and now, each part of the segment is represented by its own list of such features. Finally, the representations of all the small parts of the segment are strung together into a time series. This time series becomes the representation of that segment."},{"metadata":{"trusted":true},"cell_type":"code","source":"def transform_ts(ts, n_dim=160, min_max=(-1,1)):\n    ts_std = min_max_transfer(ts, min_value=MIN_NUM, max_value=MAX_NUM)\n    bucket_size = int(SIGNAL_LEN / n_dim)\n    new_ts = []\n    for i in range(0, SIGNAL_LEN, bucket_size):\n        ts_range = ts_std[i:i + bucket_size]\n        mean = ts_range.mean()\n        std = ts_range.std()\n        std_top = mean + std\n        std_bot = mean - std\n        percentil_calc = ts_range.quantile([0, 0.01, 0.25, 0.50, 0.75, 0.99, 1])\n        max_range = ts_range.quantile(1) - ts_range.quantile(0)\n        relative_percentile = percentil_calc - mean\n        new_ts.append(np.concatenate([np.asarray([mean, std, std_top, std_bot, max_range]), percentil_calc, relative_percentile]))\n    return np.asarray(new_ts)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Prepare the final signal features\nThe time series representations of all segments of the signal are calulated and concatenated. This results in a 3D tensor containing the time series representations of all segments in the signal."},{"metadata":{"trusted":true},"cell_type":"code","source":"def prepare_data(start, end):\n    train = pd.DataFrame(np.transpose(signals[int(start):int(end)]))\n    X = []\n    for id_measurement in tqdm(train.index[int(start):int(end)]):\n        X_signal = transform_ts(train[id_measurement])\n        X.append(X_signal)\n    X = np.asarray(X)\n    return X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Implement the feature generation process"},{"metadata":{"trusted":true},"cell_type":"code","source":"X = []\n\ndef load_all():\n    total_size = len(signals)\n    for start, end in [(0, int(total_size))]:\n        X_temp = prepare_data(start, end)\n        X.append(X_temp)\n        \nload_all()\nX = np.concatenate(X)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here is the shape of X. There are a total of 4194 segments. Each segment is divided into 161 parts and each part is represented by a list of 19 features (mean, stddev etc). Therefore, X is a 3D tensor with shape (4194, 161, 19)."},{"metadata":{"trusted":true},"cell_type":"code","source":"X.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Flattening the features and doing basic EDA with seaborn"},{"metadata":{},"cell_type":"markdown","source":"Now, we flatten the 2D tensors associated with each segment into 1D arrays. Now, each data point (segment) is represented by a 1D array.\n\nHere are some flattened 1D arrays (with sparse selection) visualized with **matplotlib**."},{"metadata":{"trusted":true},"cell_type":"code","source":"shape = X.shape\nnew_signals = X.reshape((shape[0], shape[1]*shape[2]))\n\nsparse_signals = []\nfor i in range(3):\n    sparse_signal = []\n    for j in range(len(new_signals[i])):\n        if j % 3 == 0:\n            sparse_signal.append(new_signals[i][j])\n    sparse_signals.append(sparse_signal)\n\nplt.plot(sparse_signals[0], 'mediumseagreen')\nplt.show()\nplt.plot(sparse_signals[1], 'seagreen')\nplt.show()\nplt.plot(sparse_signals[2], 'green')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Spectral Entropy\nThe spectral entropy is a method of calculating the complexity or entropy (disorderliness) of a time series. It's defined to be the Shannon Entropy of the Power Spectral Density (PSD) of the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"def spectral_entropy(x, sf, method='fft', nperseg=None, normalize=False):\n    \"\"\"Spectral Entropy.\n    Parameters\n    ----------\n    x : list or np.array\n        One-dimensional time series of shape (n_times)\n    sf : float\n        Sampling frequency\n    method : str\n        Spectral estimation method ::\n        'fft' : Fourier Transform (via scipy.signal.periodogram)\n        'welch' : Welch periodogram (via scipy.signal.welch)\n    nperseg : str or int\n        Length of each FFT segment for Welch method.\n        If None, uses scipy default of 256 samples.\n    normalize : bool\n        If True, divide by log2(psd.size) to normalize the spectral entropy\n        between 0 and 1. Otherwise, return the spectral entropy in bit.\n    Returns\n    -------\n    se : float\n        Spectral Entropy\n    Notes\n    -----\n    Spectral Entropy is defined to be the Shannon Entropy of the Power\n    Spectral Density (PSD) of the data:\n    .. math:: H(x, sf) =  -\\\\sum_{f=0}^{f_s/2} PSD(f) log_2[PSD(f)]\n    Where :math:`PSD` is the normalised PSD, and :math:`f_s` is the sampling\n    frequency.\n    References\n    ----------\n    .. [1] Inouye, T. et al. (1991). Quantification of EEG irregularity by\n       use of the entropy of the power spectrum. Electroencephalography\n       and clinical neurophysiology, 79(3), 204-210.\n    Examples\n    --------\n    1. Spectral entropy of a pure sine using FFT\n        >>> from entropy import spectral_entropy\n        >>> import numpy as np\n        >>> sf, f, dur = 100, 1, 4\n        >>> N = sf * duration # Total number of discrete samples\n        >>> t = np.arange(N) / sf # Time vector\n        >>> x = np.sin(2 * np.pi * f * t)\n        >>> print(np.round(spectral_entropy(x, sf, method='fft'), 2)\n            0.0\n    2. Spectral entropy of a random signal using Welch's method\n        >>> from entropy import spectral_entropy\n        >>> import numpy as np\n        >>> np.random.seed(42)\n        >>> x = np.random.rand(3000)\n        >>> print(spectral_entropy(x, sf=100, method='welch'))\n            9.939\n    3. Normalized spectral entropy\n        >>> print(spectral_entropy(x, sf=100, method='welch', normalize=True))\n            0.995\n    \"\"\"\n    x = np.array(x)\n    # Compute and normalize power spectrum\n    if method == 'fft':\n        _, psd = periodogram(x, sf)\n    elif method == 'welch':\n        _, psd = welch(x, sf, nperseg=nperseg)\n    psd_norm = np.divide(psd, psd.sum())\n    se = -np.multiply(psd_norm, np.log2(psd_norm)).sum()\n    if normalize:\n        se /= np.log2(psd_norm.size)\n    return se","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spectral_entropies = np.array([spectral_entropy(new_signal, sf=100, method='fft') for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=spectral_entropies, y=targets, kind='kde', color='blueviolet')\nplot.set_axis_labels('spectral_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) along a line with negative slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=spectral_entropies, y=targets, kind='hex', color='blueviolet')\nplot.set_axis_labels('spectral_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a negatively-sloped line."},{"metadata":{},"cell_type":"markdown","source":"#### Scatterplot with line of best fit"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=spectral_entropies, y=targets, kind='reg', color='blueviolet')\nplot.set_axis_labels('spectral_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear negative slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **negative correlation** between the spectral entropy of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### Sample Entropy\nSample entropy is a modification of approximate entropy (which I used in my previous kernel), generally used for assessing the complexity of physiological time-series signals. It has two advantages over approximate entropy: data length independence and a relatively trouble-free implementation. Large values indicate high complexity whereas smaller values characterize more self-similar and regular signals."},{"metadata":{"trusted":true},"cell_type":"code","source":"@jit('f8(f8[:], i4, f8)', nopython=True)\ndef _numba_sampen(x, mm=2, r=0.2):\n    \"\"\"\n    Fast evaluation of the sample entropy using Numba.\n    \"\"\"\n    n = x.size\n    n1 = n - 1\n    mm += 1\n    mm_dbld = 2 * mm\n\n    # Define threshold\n    r *= x.std()\n\n    # initialize the lists\n    run = [0] * n\n    run1 = run[:]\n    r1 = [0] * (n * mm_dbld)\n    a = [0] * mm\n    b = a[:]\n    p = a[:]\n\n    for i in range(n1):\n        nj = n1 - i\n\n        for jj in range(nj):\n            j = jj + i + 1\n            if abs(x[j] - x[i]) < r:\n                run[jj] = run1[jj] + 1\n                m1 = mm if mm < run[jj] else run[jj]\n                for m in range(m1):\n                    a[m] += 1\n                    if j < n1:\n                        b[m] += 1\n            else:\n                run[jj] = 0\n        for j in range(mm_dbld):\n            run1[j] = run[j]\n            r1[i + n * j] = run[j]\n        if nj > mm_dbld - 1:\n            for j in range(mm_dbld, nj):\n                run1[j] = run[j]\n\n    m = mm - 1\n\n    while m > 0:\n        b[m] = b[m - 1]\n        m -= 1\n\n    b[0] = n * n1 / 2\n    a = np.array([float(aa) for aa in a])\n    b = np.array([float(bb) for bb in b])\n    p = np.true_divide(a, b)\n    return -log(p[-1])\n\ndef sample_entropy(x, order=2, metric='chebyshev'):\n    \"\"\"Sample Entropy.\n    Parameters\n    ----------\n    x : list or np.array\n        One-dimensional time series of shape (n_times)\n    order : int (default: 2)\n        Embedding dimension.\n    metric : str (default: chebyshev)\n        Name of the metric function used with KDTree. The list of available\n        metric functions is given by: `KDTree.valid_metrics`.\n    Returns\n    -------\n    se : float\n        Sample Entropy.\n    Notes\n    -----\n    Sample entropy is a modification of approximate entropy, used for assessing\n    the complexity of physiological time-series signals. It has two advantages\n    over approximate entropy: data length independence and a relatively\n    trouble-free implementation. Large values indicate high complexity whereas\n    smaller values characterize more self-similar and regular signals.\n    Sample entropy of a signal :math:`x` is defined as:\n    .. math:: H(x, m, r) = -log\\\\frac{C(m + 1, r)}{C(m, r)}\n    where :math:`m` is the embedding dimension (= order), :math:`r` is\n    the radius of the neighbourhood (default = :math:`0.2 * \\\\text{std}(x)`),\n    :math:`C(m + 1, r)` is the number of embedded vectors of length\n    :math:`m + 1` having a Chebyshev distance inferior to :math:`r` and\n    :math:`C(m, r)` is the number of embedded vectors of length\n    :math:`m` having a Chebyshev distance inferior to :math:`r`.\n    Note that if metric == 'chebyshev' and x.size < 5000 points, then the\n    sample entropy is computed using a fast custom Numba script. For other\n    metric types or longer time-series, the sample entropy is computed using\n    a code from the mne-features package by Jean-Baptiste Schiratti\n    and Alexandre Gramfort (requires sklearn).\n    References\n    ----------\n    .. [1] Richman, J. S. et al. (2000). Physiological time-series analysis\n           using approximate entropy and sample entropy. American Journal of\n           Physiology-Heart and Circulatory Physiology, 278(6), H2039-H2049.\n    Examples\n    --------\n    1. Sample entropy with order 2.\n        >>> from entropy import sample_entropy\n        >>> import numpy as np\n        >>> np.random.seed(1234567)\n        >>> x = np.random.rand(3000)\n        >>> print(sample_entropy(x, order=2))\n            2.192\n    2. Sample entropy with order 3 using the Euclidean distance.\n        >>> from entropy import sample_entropy\n        >>> import numpy as np\n        >>> np.random.seed(1234567)\n        >>> x = np.random.rand(3000)\n        >>> print(sample_entropy(x, order=3, metric='euclidean'))\n            2.725\n    \"\"\"\n    x = np.asarray(x, dtype=np.float64)\n    if metric == 'chebyshev' and x.size < 5000:\n        return _numba_sampen(x, mm=order, r=0.2)\n    else:\n        phi = _app_samp_entropy(x, order=order, metric=metric,\n                                approximate=False)\n        return -np.log(np.divide(phi[1], phi[0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_entropies = np.array([sample_entropy(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=sample_entropies, y=targets, kind='kde', color='mediumvioletred')\nplot.set_axis_labels('sample_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) around a line with negative slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=sample_entropies, y=targets, kind='hex', color='mediumvioletred')\nplot.set_axis_labels('sample_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a negatively-sloped line."},{"metadata":{},"cell_type":"markdown","source":"#### Scatterplot with line of best fit"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=sample_entropies, y=targets, kind='reg', color='mediumvioletred')\nplot.set_axis_labels('sample_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear negative slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **negative correlation** between the sample entropy of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### Detrended Fluctuation\nDetrended fluctuation analysis (DFA) is used to find long-term statistical dependencies in time series. It is another excellent way to measure the entropy or complexity of a signal. For more details, please refer to the excellent documentation of the nold Python package by Christopher Scholzel, from which this function is taken: https://cschoel.github.io/nolds/nolds.html#detrended-fluctuation-analysis."},{"metadata":{"trusted":true},"cell_type":"code","source":"@jit('UniTuple(float64, 2)(float64[:], float64[:])', nopython=True)\ndef _linear_regression(x, y):\n    \"\"\"Fast linear regression using Numba.\n    Parameters\n    ----------\n    x, y : ndarray, shape (n_times,)\n        Variables\n    Returns\n    -------\n    slope : float\n        Slope of 1D least-square regression.\n    intercept : float\n        Intercept\n    \"\"\"\n    n_times = x.size\n    sx2 = 0\n    sx = 0\n    sy = 0\n    sxy = 0\n    for j in range(n_times):\n        sx2 += x[j] ** 2\n        sx += x[j]\n        sxy += x[j] * y[j]\n        sy += y[j]\n    den = n_times * sx2 - (sx ** 2)\n    num = n_times * sxy - sx * sy\n    slope = num / den\n    intercept = np.mean(y) - slope * np.mean(x)\n    return slope, intercept\n\n\n@jit('i8[:](f8, f8, f8)', nopython=True)\ndef _log_n(min_n, max_n, factor):\n    \"\"\"\n    Creates a list of integer values by successively multiplying a minimum\n    value min_n by a factor > 1 until a maximum value max_n is reached.\n    Used for detrended fluctuation analysis (DFA).\n    Function taken from the nolds python package\n    (https://github.com/CSchoel/nolds) by Christopher Scholzel.\n    Parameters\n    ----------\n    min_n (float):\n        minimum value (must be < max_n)\n    max_n (float):\n        maximum value (must be > min_n)\n    factor (float):\n       factor used to increase min_n (must be > 1)\n    Returns\n    -------\n    list of integers:\n        min_n, min_n * factor, min_n * factor^2, ... min_n * factor^i < max_n\n        without duplicates\n    \"\"\"\n    max_i = int(floor(log(1.0 * max_n / min_n) / log(factor)))\n    ns = [min_n]\n    for i in range(max_i + 1):\n        n = int(floor(min_n * (factor ** i)))\n        if n > ns[-1]:\n            ns.append(n)\n    return np.array(ns, dtype=np.int64)\n\n@jit('f8(f8[:])', nopython=True)\ndef _dfa(x):\n    \"\"\"\n    Utility function for detrended fluctuation analysis\n    \"\"\"\n    N = len(x)\n    nvals = _log_n(4, 0.1 * N, 1.2)\n    walk = np.cumsum(x - x.mean())\n    fluctuations = np.zeros(len(nvals))\n\n    for i_n, n in enumerate(nvals):\n        d = np.reshape(walk[:N - (N % n)], (N // n, n))\n        ran_n = np.array([float(na) for na in range(n)])\n        d_len = len(d)\n        slope = np.empty(d_len)\n        intercept = np.empty(d_len)\n        trend = np.empty((d_len, ran_n.size))\n        for i in range(d_len):\n            slope[i], intercept[i] = _linear_regression(ran_n, d[i])\n            y = np.zeros_like(ran_n)\n            # Equivalent to np.polyval function\n            for p in [slope[i], intercept[i]]:\n                y = y * ran_n + p\n            trend[i, :] = y\n        # calculate standard deviation (fluctuation) of walks in d around trend\n        flucs = np.sqrt(np.sum((d - trend) ** 2, axis=1) / n)\n        # calculate mean fluctuation over all subsequences\n        fluctuations[i_n] = flucs.sum() / flucs.size\n\n    # Filter zero\n    nonzero = np.nonzero(fluctuations)[0]\n    fluctuations = fluctuations[nonzero]\n    nvals = nvals[nonzero]\n    if len(fluctuations) == 0:\n        # all fluctuations are zero => we cannot fit a line\n        dfa = np.nan\n    else:\n        dfa, _ = _linear_regression(np.log(nvals), np.log(fluctuations))\n    return dfa\n\n\ndef detrended_fluctuation(x):\n    \"\"\"\n    Detrended fluctuation analysis (DFA).\n    Parameters\n    ----------\n    x : list or np.array\n        One-dimensional time-series.\n    Returns\n    -------\n    dfa : float\n        the estimate alpha for the Hurst parameter:\n        alpha < 1: stationary process similar to fractional Gaussian noise\n        with H = alpha\n        alpha > 1: non-stationary process similar to fractional Brownian\n        motion with H = alpha - 1\n    Notes\n    -----\n    Detrended fluctuation analysis (DFA) is used to find long-term statistical\n    dependencies in time series.\n    The idea behind DFA originates from the definition of self-affine\n    processes. A process :math:`X` is said to be self-affine if the standard\n    deviation of the values within a window of length n changes with the window\n    length factor L in a power law:\n    .. math:: \\\\text{std}(X, L * n) = L^H * \\\\text{std}(X, n)\n    where :math:`\\\\text{std}(X, k)` is the standard deviation of the process\n    :math:`X` calculated over windows of size :math:`k`. In this equation,\n    :math:`H` is called the Hurst parameter, which behaves indeed very similar\n    to the Hurst exponent.\n    For more details, please refer to the excellent documentation of the nolds\n    Python package by Christopher Scholzel, from which this function is taken:\n    https://cschoel.github.io/nolds/nolds.html#detrended-fluctuation-analysis\n    Note that the default subseries size is set to\n    entropy.utils._log_n(4, 0.1 * len(x), 1.2)). The current implementation\n    does not allow to manually specify the subseries size or use overlapping\n    windows.\n    The code is a faster (Numba) adaptation of the original code by Christopher\n    Scholzel.\n    References\n    ----------\n    .. [1] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons,\n           H. E. Stanley, and A. L. Goldberger, “Mosaic organization of\n           DNA nucleotides,” Physical Review E, vol. 49, no. 2, 1994.\n    .. [2] R. Hardstone, S.-S. Poil, G. Schiavone, R. Jansen,\n           V. V. Nikulin, H. D. Mansvelder, and K. Linkenkaer-Hansen,\n           “Detrended fluctuation analysis: A scale-free view on neuronal\n           oscillations,” Frontiers in Physiology, vol. 30, 2012.\n    Examples\n    --------\n        >>> import numpy as np\n        >>> from entropy import detrended_fluctuation\n        >>> np.random.seed(123)\n        >>> x = np.random.rand(100)\n        >>> print(detrended_fluctuation(x))\n            0.761647725305623\n    \"\"\"\n    x = np.asarray(x, dtype=np.float64)\n    return _dfa(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"detrended_fluctuations = np.array([detrended_fluctuation(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=detrended_fluctuations, y=targets, kind='kde', color='mediumblue')\nplot.set_axis_labels('detrended_fluctuation', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) around a line with negative slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=detrended_fluctuations, y=targets, kind='hex', color='mediumblue')\nplot.set_axis_labels('detrended_fluctuation', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a negatively-sloped line."},{"metadata":{},"cell_type":"markdown","source":"#### Scatterplot with line of best fit"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot = sns.jointplot(x=detrended_fluctuations, y=targets, kind='reg', color='mediumblue')\nplot.set_axis_labels('detrended_fluctuation', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear negative slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **negative correlation** between the detrended fluctuation of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"## Conclusion\nEverywhere, we see that the features measuring the complexity or roughess of the curve (entropy, fractal dimesion, detrended fluctuation etc) are negatively correlated with the time left for the earthquake to occur. **This makes sense, because as the earthquake gets closer, and the time to failure decreases, the seismic signal should become more complex, rough and unpredictable**. This is down to the fact that the seismc signal's energy and amplitude generally increase rapidly as the earthquake comes closer. Therefore, these features along with energy-based features can be useful in predicting the time left for an earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### That's it ! Thanks for reading this kernel :)"},{"metadata":{},"cell_type":"markdown","source":"### Please post your feedback and suggestions in the comments."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.4"}},"nbformat":4,"nbformat_minor":1}