{
 "metadata": {
  "name": "",
  "signature": "sha256:c6b1c463d94ec780d4ea419d3046a5fd4cab4a33c35f82ef3be3bf60cfe6e1dd"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Scikit-learn's batch logistic regression doesn't support a binomial response -- where your data is (successes, trials) rather than (success or failure). This can be a drag when you have a lot of trials with the exact same features, as your design matrix will be repetitive and much larger than necessary.\n",
      "\n",
      "It isn't too bad to do it yourself though, which I'll show here. But first, my notation:\n",
      "\n",
      "* $n$ - number of trials\n",
      "* $k$ - number of successes\n",
      "* $X$ - design matrix\n",
      "* $\\lambda$ - regularization parameter; I use `lambd` in the code because `lambda` is reserved in python\n",
      "\n",
      "In the code, I treat the intercept separately and don't regularize it, though it's pretty straightforward to modify things and include it as a normal feature."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.optimize import fmin_l_bfgs_b\n",
      "from numpy import log, exp, hstack\n",
      "\n",
      "def logistic(x):\n",
      "    return 1/(1+exp(-x))\n",
      " \n",
      "def log_loss(y_orig, n, k, eps=1e-15):\n",
      "    y = y_orig.copy()\n",
      "    y[y < eps] = eps\n",
      "    y[y > 1 - eps] = 1 - eps\n",
      "    x = -k * log(y) - (n - k) * log(1 - y)\n",
      "    return x.sum()/n.sum()\n",
      " \n",
      "def lbfgs_func(x, X, n, k, lambd=1.0):\n",
      "    intercept, theta = x[0], x[1:]\n",
      "    y = logistic(X.dot(theta) + intercept)\n",
      "    penalty = lambd/2*(theta*theta).sum()/n.sum()\n",
      "    obj = log_loss(y, n, k) + penalty\n",
      " \n",
      "    r = n*logistic(X.dot(theta)+intercept) - k\n",
      "    grad_intercept = r.sum()/n.sum()\n",
      "    grad_coefs = (X.T.dot(r) + lambd*theta)/n.sum()\n",
      "    grad = hstack([grad_intercept, grad_coefs])\n",
      " \n",
      "    return obj, grad"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To use the code, you need your design matrix $X$ and a vector each of trials and successes, $n$ and $k$. It'll look something like this:\n",
      "\n",
      "```python\n",
      "x0 = np.zeros(X.shape[1]+1)\n",
      "x_optimal, objective, info = fmin_l_bfgs_b(lbfgs_func, x0, args=(X, n, k))\n",
      "intercept, theta = x_optimal[0], x_optimal[1:]\n",
      "```\n",
      "\n",
      "Now `intercept` and `theta` will have your optimal parameters. To make predictions on new data, do something like this:\n",
      "\n",
      "```python\n",
      "y_pred = logistic(X_new.dot(theta) + intercept)\n",
      "```\n",
      "\n",
      "Note that for large problems you might want to play with the defaults for fmin_l_bfgs_b; in particular I've found changing the factr keyword argument to be helpful."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Example\n",
      "\n",
      "Two features with four training samples each. The first feature got 3/4 successes, the second only one out of four. You can see that scikit-learn's answer agrees with mine."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from sklearn import linear_model"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.matrix('1 0; 1 0; 1 0; 1 0; 0 1; 0 1; 0 1; 0 1').A\n",
      "y = np.array([1, 1, 1, 0, 1, 0, 0, 0])\n",
      "\n",
      "C = 10\n",
      "clf = linear_model.LogisticRegression(C=C)\n",
      "clf.fit(X, y)\n",
      "clf.intercept_, clf.coef_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "(array([ 0.]), array([[ 0.97281217, -0.97281217]]))"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.array([[1, 0], [0, 1]])\n",
      "n = np.array([4, 4])\n",
      "k = np.array([3, 1])\n",
      "\n",
      "lambd = 1/C\n",
      "x0 = np.zeros(X.shape[1]+1)\n",
      "x_optimal, *_ = fmin_l_bfgs_b(lbfgs_func, x0, args=(X, n, k, lambd))\n",
      "intercept, theta = x_optimal[0], x_optimal[1:]\n",
      "intercept, theta"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "(-5.1209353873023446e-16, array([ 0.97280393, -0.97280393]))"
       ]
      }
     ],
     "prompt_number": 4
    }
   ],
   "metadata": {}
  }
 ]
}