Analisys_Code_Maximum_Likelihood.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "F0OiUYETtNBU"
      },
      "outputs": [],
      "source": [
        "# VALÉRIO DE ARRUDA, M. (marcel.a@ufabc.edu.br)\n",
        "\n",
        "#======================================================\n",
        "#===========                    =======================\n",
        "#=========== MAXIMUM LIKELIHOOD =======================\n",
        "#===========                    =======================\n",
        "#======================================================\n",
        "\n",
        "# Directory indication\n",
        "#import os\n",
        "#os.getcwd()\n",
        "#os.chdir(\"/Users/your_username\")\n",
        "\n",
        "# Libraries\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot  as plt\n",
        "import datetime\n",
        "from datetime import timedelta, timezone, datetime\n",
        "from scipy import stats\n",
        "from scipy.optimize import minimize\n",
        "from scipy.optimize import curve_fit\n",
        "import pylab as py\n",
        "import matplotlib.patches as mpatches\n",
        "\n",
        "# Data reading (.csv)\n",
        "name = '/content/S_124.csv' # OpenSesame/JATOS Table (.csv) path\n",
        "df = pd.read_csv(name, header = 0)\n",
        "df\n",
        "# Document/participant ID\n",
        "name_file = (name.split(\".\")[0])\n",
        "name_file = name_file.replace(\"/content/\", \"\")\n",
        "name_file\n",
        "\n",
        "# Data pre-processing\n",
        "post_practice = df.response_resp_log.size - 70\n",
        "error_practice = post_practice - 20\n",
        "perc_answer = pd.crosstab(df.response_resp_log[post_practice:], df.tempo_log, margins=True, normalize = 'columns')\n",
        "time_stm = np.array([500.0, 630.0, 793.5, 1000.0, 1260.0, 1587.0, 2000])\n",
        "perc_long = perc_answer[0:1].values\n",
        "perc_long = np.delete(perc_long,[7])\n",
        "\n",
        "# SIGMOID CURVE ADJUSTMENT WITH MAXIMUM LIKELIHOOD =======================================================\n",
        "def sigmoid(params):\n",
        "    k = params[0]    # parameter a from the previous code\n",
        "    x0 = params[1]   \n",
        "    sd = params[2]   # parameter b from the previous code\n",
        "\n",
        "    yPred = 1 / (1+ np.exp(-k*(time_stm-x0)))\n",
        "            \n",
        "    # Calculate negative log likelihood\n",
        "    LL = -np.sum(stats.norm.logpdf(perc_long, loc=yPred, scale=sd))\n",
        "\n",
        "    return(LL)\n",
        "\n",
        "# Definition of initial parameters\n",
        "# initParams = [k, inflection point, sd]\n",
        "# I'm still trying to figure out what k and sd are.\n",
        "\n",
        "# sd according to the reference code is the Standard Deviation.\n",
        "# but here in practice it seems to be related to the slope of the fitted curve.\n",
        "\n",
        "# here are the guesses for the initial parameters.\n",
        "# the second parameter is a guess for the inflection point (I have used geometric or arithmetic mean)\n",
        "initParams = [0.0, 1000, 0.005]\n",
        "\n",
        "# minimization of the distance between the initial parameters and the parameters of the sigmoid obtained\n",
        "results = minimize(sigmoid, initParams, method='Nelder-Mead')\n",
        "print(\"Parametros Iniciais [k, ponto de inflexão, sd]:\\n\",results.x)\n",
        "estParms = results.x\n",
        "\n",
        "rng = np.arange(time_stm[0],time_stm[-1],0.1)\n",
        "\n",
        "yOut = yPred = 1 / (1+ np.exp(-estParms[0]*(rng-estParms[1])))\n",
        "\n",
        "\n",
        "# Bisection Point===========================================\n",
        "\n",
        "pnt_bisec = round(estParms[1],4)\n",
        "\n",
        "# Differential Threshold====================================\n",
        "\n",
        "    #0,25\n",
        "T1_ = estParms[1] - (np.log(3)/estParms[0])\n",
        "\n",
        "    #0,75\n",
        "T2_ = estParms[1] + (np.log(3)/estParms[0])\n",
        "\n",
        "lim_dif_ = round(((T2_-T1_)/2),4)\n",
        "\n",
        "# Weber's Ratio==============================================\n",
        "\n",
        "WR = round(((T2_-T1_)/pnt_bisec),4)\n",
        "\n",
        "# Print Results==============================================\n",
        "r = {'Participant': [name_file], 'BP': [pnt_bisec], 'LD': [lim_dif_], 'WR':[WR]}\n",
        "Results = pd.DataFrame(data=r)\n",
        "print(\"\\nResultados:\\n\")\n",
        "display(Results)\n",
        "\n",
        "# GRAPH COMMANDs=============================================\n",
        "print(\"\\nPlot:\\n\")\n",
        "\n",
        "# Import seaborn\n",
        "import seaborn as sns\n",
        "\n",
        "# Apply the default theme\n",
        "sns.set_theme()\n",
        "sns.set_style(\"ticks\")\n",
        "\n",
        "\n",
        "rng = np.arange(time_stm[0],time_stm[-1],0.1)\n",
        "pnt_bisec = round(estParms[1],4)\n",
        "\n",
        "plt.figure(figsize = (4,3),dpi = 200)\n",
        "plt.plot(rng, yOut, label=\"Curve-fit [Max Likelihood]\", linestyle ='-',lw = 3,color='k')\n",
        "plt.plot(time_stm, perc_long, color='k', marker='o',linestyle = '',ms = 7,markerfacecolor = 'w')\n",
        "plt.ylabel('P(Long)',fontsize = 10)\n",
        "plt.xlabel('Duration (ms)',fontsize = 10)\n",
        "plt.ylim(-0.1,1.1)\n",
        "plt.yticks(np.arange(0, 1.1, step=0.25))\n",
        "plt.grid()\n",
        "\n",
        "#plt.savefig('BISSEC_{0}'.format(name_file), transparent=True,bbox_inches='tight',format='png')\n",
        "plt.show()"
      ]
    }
  ]
}