TransferEntropy_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import random\n",
    "import time\n",
    "import math\n",
    "import itertools\n",
    "\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def joint_probability_new_version(k,l,h, a, b, lbl_a, lbl_b):\n",
    "    '''\n",
    "    k B time horizon\n",
    "    l A time horizon\n",
    "    h instant in the future of serie B\n",
    "\n",
    "    a, b array type\n",
    "\n",
    "    lbl_a - label of a series\n",
    "    lbl_b - label of b series\n",
    "    '''\n",
    "        \n",
    "    \n",
    "    numStates=2**(k+l+1)\n",
    "    combinations = list(map(list, itertools.product([0, 1], repeat=k+l+1)))\n",
    "    prob_cnjt = np.zeros(numStates)\n",
    "\n",
    "    matrix_nova = np.matrix([b[1:],b[:-1],a[:-1]]).T\n",
    "    df = pd.DataFrame(matrix_nova, columns = ['b_ftr', lbl_b, lbl_a])\n",
    "    gpd = df.groupby(['b_ftr', lbl_b, lbl_a], as_index=False).size().reset_index(name='Count')\n",
    "    total = sum(gpd['Count'])\n",
    "    \n",
    "    for i in np.arange(0,gpd.shape[0]):\n",
    "        comb = [e for e in gpd.iloc[i][0:3].values.tolist()]\n",
    "        idx = combinations.index(comb)\n",
    "        prob_cnjt[idx] = gpd.iloc[i]['Count']/total\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def joint_probability(k,l,h, a, b):\n",
    "    '''\n",
    "        k B time horizon\n",
    "        l A time horizon\n",
    "        h instant in the future of serie B\n",
    "        \n",
    "        a, b array type'''\n",
    "\n",
    "    #Alarm Series A (cause), B (effect), same len\n",
    "    #tested\n",
    "    sizeSeries = a.size\n",
    "    transEntropy = 0\n",
    "    numStates = 2**(k + l  + 1)\n",
    "    combinations = list(map(list, itertools.product([0, 1], repeat=k+l+1)))\n",
    "    counting = np.zeros(numStates)\n",
    "    prob_cnjt = np.zeros(numStates)\n",
    "    a_prob_ind = []\n",
    "    b_prob_ind = []\n",
    "    #joitn probability p(i_sub_t+1), i_sub_t**k, j_sub_t**l)\n",
    "    inicio = np.max([k,l]) - 1\n",
    "    for i in np.arange(inicio, sizeSeries - h):\n",
    "        for hk in np.arange(0,k):\n",
    "                b_prob_ind.append(b[i - hk])\n",
    "        for hl in np.arange(0,l):\n",
    "                a_prob_ind.append(a[i - hl])\n",
    "\n",
    "        #print(a.size, b.size, a.size -1)     \n",
    "        ab = [b[i + h]] + b_prob_ind + a_prob_ind \n",
    "        index_comb = combinations.index(ab)\n",
    "        counting[index_comb] = counting[index_comb] + 1\n",
    "\n",
    "        a_prob_ind = []\n",
    "        b_prob_ind = []\n",
    "\n",
    "    total = sum(counting)\n",
    "  \n",
    "    prob_cnjt = counting/total\n",
    "     \n",
    "    return prob_cnjt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Joint probability evaluation p(i_t+h, i_t**k)\n",
    "#tested\n",
    "def joint_prob_ih_ik(k,l, joint_prob_ih_ik_jl):\n",
    "    states_ith_ik = list(map(list, itertools.product([0, 1], repeat=k + 1)))\n",
    "    combinations = list(map(list, itertools.product([0, 1], repeat=k+l+1))) \n",
    "    p_jnt_ith_ik = np.zeros(2**(k+1))\n",
    "    \n",
    "    for i, state in enumerate(states_ith_ik):\n",
    "        for j, comb in enumerate(combinations):\n",
    "            if comb[0:k+1] == state:\n",
    "                p_jnt_ith_ik[i] = p_jnt_ith_ik[i] + joint_prob_ih_ik_jl[j]\n",
    "    return p_jnt_ith_ik"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conditional_prob(k,l,joint_prob):\n",
    "    states = list(map(list, itertools.product([0, 1], repeat=k+l)))\n",
    "    combinations = list(map(list, itertools.product([0, 1], repeat=k+l+1)))\n",
    "\n",
    "    size = int(joint_prob.size/2)\n",
    "    conditional = np.zeros(2**(k+l+1))\n",
    "\n",
    "    for i,state in enumerate(states):\n",
    "        index_zero = combinations.index([0] + state)\n",
    "        prob_zero = joint_prob[index_zero]\n",
    "\n",
    "        index_one = combinations.index([1] + state)\n",
    "        prob_one = joint_prob[index_one]\n",
    "\n",
    "        if(prob_zero + prob_one != 0):\n",
    "            conditional[i] = prob_zero/(prob_zero+ prob_one)\n",
    "            conditional[i + 2**(k+l)] = prob_one/(prob_zero+ prob_one)\n",
    "    return conditional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Division of the conditionals in log2 \n",
    "#tested\n",
    "def conditional_div(k,l,conditional_num, conditional_den):\n",
    "    combinations = list(map(list, itertools.product([0, 1], repeat=k+l+1)))\n",
    "    conditional_division = np.zeros(conditional_num.size)\n",
    "    states_den = list(map(list, itertools.product([0, 1], repeat=1+k)))\n",
    "    for j, comb in enumerate(combinations):\n",
    "        if(conditional_den[states_den.index(comb[0:k+1])] != 0):\n",
    "            conditional_division[j] = conditional_num[j]/conditional_den[states_den.index(comb[0:k+1])]            \n",
    "    return conditional_division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transfer entropy final evaluation\n",
    "def te(k,l,h,a,b):\n",
    "    '''\n",
    "        transentropy a->b\n",
    "        te(k,l,h,a,b)\n",
    "        k - dimension of b, number of samples of the past of b\n",
    "        l - dimension of a, number of samples of the passt of a\n",
    "        h -> instant in the future of b\n",
    "    '''\n",
    "    joint_p_ih_ik_jl = joint_probability(k,l,h,a,b)\n",
    "    \n",
    "    joint_p_ih_ik = joint_prob_ih_ik(k,l, joint_p_ih_ik_jl)\n",
    "    conditional_num = conditional_prob(k,l,joint_p_ih_ik_jl)\n",
    "    conditional_den = conditional_prob(k,0, joint_p_ih_ik)    \n",
    "    div = conditional_div(k,l,conditional_num, conditional_den)\n",
    "    \n",
    "    #log2 from the division of the conditionals -> #p(i_sub_t+h|i_sub_t**k, j_sub_t**l) /p(i_sub_t+h|i_t**k)\n",
    "    \n",
    "    log2_div_cond = np.log2(div[div!=0])\n",
    "    te = np.sum(joint_p_ih_ik_jl[div!=0]*log2_div_cond)\n",
    "    return te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transferEntropy_case(df, k, l, h):\n",
    "    \n",
    "    '''Evaluate Transfer entropy for a dataframe of variables'''\n",
    "    start = time.clock()\n",
    "    transEntropy = np.zeros([df.columns.size,df.columns.size])\n",
    "    sigValues =  np.zeros([df.columns.size,df.columns.size])\n",
    "    for i in np.arange(0, df.columns.size):\n",
    "        for j in np.arange(0, df.columns.size):\n",
    "            print('trans ', df.columns[i], df.columns[j])\n",
    "            if(j != i + df.columns.size/2 and j!=i and j != i - df.columns.size/2):\n",
    "                transEntropy[i][j] = te(k,l,h,df[df.columns[i]],\n",
    "                                         df[df.columns[j]])\n",
    "            clear_output()\n",
    "    end = time.clock()   \n",
    "    \n",
    "    print(end - start)\n",
    "    return transEntropy  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#joint probablity for functions test\n",
    "#joint_p_ih_ik_jl = np.array([0.97322404,0.00546448,0.00491803,0,0,0.00546448, 0.00546448, 0.00546448])\n",
    "\n",
    "#aproximate results for this test\n",
    "\n",
    "#p(ith, ik)\n",
    "#jnt_p_ih_ik = [0.97868852,0.00491803,0.00546448,0.0109286] \n",
    "\n",
    "#p(i_t+h|i**k, j**l)\n",
    "#cond_p_ih_ik_jl =  [1,0.5,0.4736841094123,0,0,0.5, 0.52631589085076,1]\n",
    "\n",
    "#p(i_th|i_k)\n",
    "#cond_p_ih_ik = [0.994711793480152,0.31035179088550,0.0552469991962,0.68964820911449]"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}