likelihood.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sparse\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.sparse import lil_matrix\n",
    "from scipy.sparse import csc_matrix\n",
    "from collections import defaultdict\n",
    "from sklearn.preprocessing import normalize\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DATA AS MATRICES\n",
    "# S = csr_matrix(np.array([[0, 1, 1, 0],\n",
    "#                          [1, 0, 0, 1],\n",
    "#                          [1, 0, 0, 0],\n",
    "#                          [0, 1, 0, 0]]))\n",
    "\n",
    "# T = csr_matrix(np.array([[0, 10,10, 0],\n",
    "#                          [10, 0, 0,10],\n",
    "#                          [10, 0, 0, 0],\n",
    "#                          [0, 10, 0, 0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DATA AS GRAPHS\n",
    "# STRUCTURE\n",
    "G=nx.DiGraph()\n",
    "G.add_edge('a','y')\n",
    "G.add_edge('y','a')\n",
    "G.add_edge('a','m')\n",
    "G.add_edge('m','a')\n",
    "G.add_edge('y','y')\n",
    "\n",
    "### TRANSITIONS\n",
    "H=nx.DiGraph()\n",
    "H.add_edge('a','y',weight=10)\n",
    "H.add_edge('a','m',weight=5)\n",
    "H.add_edge('m','y',weight=1)\n",
    "\n",
    "### TO UNDIRECTED\n",
    "G = G.to_directed()\n",
    "tmp = nx.Graph()\n",
    "tmp.add_edges_from(H.edges(), weight=0)\n",
    "for u, v, d in H.edges(data=True):\n",
    "    tmp[u][v]['weight'] += d['weight']\n",
    "H = tmp\n",
    "\n",
    "### ADJACENCY MATRICES\n",
    "S = nx.adjacency_matrix(G,['y','a','m'])\n",
    "T = nx.adjacency_matrix(H,['y','a','m'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "# AIC\n",
    "##############################################################################\n",
    "def AIC(loglikelihood, params):\n",
    "    return (-2 * loglikelihood) + (2 * params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "# RANDOM WALK\n",
    "##############################################################################\n",
    "\n",
    "def random_walk_likelihood(structure, transitions, alpha=0.85):\n",
    "    params = []\n",
    "    if alpha is not None and (alpha < 0 or alpha > 1):\n",
    "        raise ValueError(\"damping factor must be between 0 and 1 (inclusive).\")\n",
    "        return None, None\n",
    "\n",
    "    params = [alpha]\n",
    "    N = structure.shape[0]\n",
    "    ### rows and columns that matter (where there is a link and transitions)\n",
    "    nnz = structure.multiply(transitions)\n",
    "    nnz = sparse.find(nnz)\n",
    "    ### probability of following sructure: alpha / outlinks\n",
    "    if alpha > 0:\n",
    "        P = csr_matrix((alpha * structure) / structure.sum(axis=1))\n",
    "        ### only cells that matter\n",
    "        P = P[nnz[0],nnz[1]]\n",
    "        ### adding teleportation\n",
    "        P += (1-alpha) / N\n",
    "    else:\n",
    "        ### adding teleportation\n",
    "        P = (1-alpha) * structure[nnz[0],nnz[1]] / N\n",
    "    ### log-likelihood\n",
    "    l = (transitions[nnz[0],nnz[1]].A1 * np.log(P.A1)).sum()\n",
    "    return l,params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "# HOP RANK\n",
    "##############################################################################\n",
    "\n",
    "def hop_counts(structure):\n",
    "    hops = csr_matrix(structure.shape, dtype=np.int16)\n",
    "    previous_hops = None\n",
    "    \n",
    "    hop = 1\n",
    "    structure.setdiag(0)\n",
    "    m = structure.copy()\n",
    "\n",
    "    for r,row in enumerate(structure):\n",
    "        hops[r,row.indices] = hop\n",
    "        \n",
    "    while m.sum() > 0:\n",
    "        \n",
    "        if previous_hops is None:\n",
    "            previous_hops = m.copy()\n",
    "        else:\n",
    "            previous_hops += m\n",
    "    \n",
    "        m = lil_matrix(m.dot(structure))\n",
    "        m.setdiag(0)\n",
    "        m -= previous_hops\n",
    "        m = (m>0).astype(np.int8)\n",
    "        m.eliminate_zeros()\n",
    "        \n",
    "        hop += 1\n",
    "        for r,row in enumerate(m):\n",
    "            hops[r,row.indices] = hop\n",
    "\n",
    "    hops.setdiag(0) \n",
    "    return hops\n",
    "\n",
    "def hop_alphas(transitions, hops):    \n",
    "    maxhops = int(hops.max())\n",
    "    total = transitions.sum()\n",
    "    alphas = np.zeros(maxhops)\n",
    "    counts = np.zeros(maxhops)\n",
    "    \n",
    "    for hop in range(maxhops):\n",
    "        transitions_in_hop = np.any(hops==hop+1).astype(np.int8).multiply(transitions)\n",
    "        counts[hop] = transitions_in_hop.sum()\n",
    "        alphas[hop] = transitions_in_hop.sum() / total\n",
    "        \n",
    "    return alphas, counts\n",
    "\n",
    "def hop_rank_likelihood(structure, transitions):\n",
    "    hops = hop_counts(structure)\n",
    "    alphas, counts = hop_alphas(transitions, hops)\n",
    "    l = 1.0\n",
    "    for r, row in enumerate(transitions):\n",
    "        if row.sum() > 0:\n",
    "            lrow = csr_matrix(np.take(alphas,hops[r].data-1)[row.indices],dtype=np.float)\n",
    "            l += (row.data * np.log(lrow.data)).sum()\n",
    "    return l,alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "# MARKOV CHAIN\n",
    "##############################################################################\n",
    "\n",
    "def markov_chain_likelihood(structure, transitions):\n",
    "    P = normalize(transitions,'l1',axis=1)\n",
    "    l = (transitions.data * np.log(P.data)).sum()\n",
    "    params = sparse.find(P.tocsr())[2].tolist()\n",
    "    return l,params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "random walk [0.0]: \n",
      "- likelihood: -32.9583686600433\n",
      "- AIC: 67.9167373200866\n",
      "\n",
      "random walk [1.0]: \n",
      "- likelihood: -17.328679513998633\n",
      "- AIC: 36.657359027997266\n",
      "\n",
      "random walk [0.85]: \n",
      "- likelihood: -19.13781445197653\n",
      "- AIC: 40.27562890395306\n",
      "\n",
      "hop rank [ 0.9375  0.0625]: \n",
      "- likelihood: -6.481333078606697\n",
      "- AIC: 16.962666157213395\n",
      "\n",
      "markov chain [0.6666666666666666, 0.16666666666666666, 0.9090909090909091, 0.8333333333333334, 0.09090909090909091, 0.3333333333333333]: \n",
      "- likelihood: -15.60207684846164\n",
      "- AIC: 43.204153696923285\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/site-packages/scipy/sparse/compressed.py:742: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "### LOG-LIKELIHOODS\n",
    "alpha = 0.\n",
    "l,params = random_walk_likelihood(S,T,alpha)\n",
    "aic = AIC(l,len(params))\n",
    "print('random walk {}: \\n- likelihood: {}\\n- AIC: {}\\n'.format(params,l,aic))\n",
    "\n",
    "alpha = 1.\n",
    "l,params = random_walk_likelihood(S,T,alpha)\n",
    "aic = AIC(l,len(params))\n",
    "print('random walk {}: \\n- likelihood: {}\\n- AIC: {}\\n'.format(params,l,aic))\n",
    "\n",
    "alpha = 0.85\n",
    "l,params = random_walk_likelihood(S,T,alpha)\n",
    "aic = AIC(l,len(params))\n",
    "print('random walk {}: \\n- likelihood: {}\\n- AIC: {}\\n'.format(params,l,aic))\n",
    "\n",
    "l,params = hop_rank_likelihood(S,T)\n",
    "aic = AIC(l,len(params))\n",
    "print('hop rank {}: \\n- likelihood: {}\\n- AIC: {}\\n'.format(params,l,aic))\n",
    "\n",
    "l,params = markov_chain_likelihood(S,T)\n",
    "aic = AIC(l,len(params))\n",
    "print('markov chain {}: \\n- likelihood: {}\\n- AIC: {}\\n'.format(params,l,aic))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}