Why is this implementation of the conditional logit gradient failing?
<p>I've written a very simple implementation of the likelihood/gradient for the conditional logit model (explained <a href="http://data.princeton.edu/wws509/notes/c6s3.html" rel="nofollow">here</a>) - the likelihood works fine, but the gradient isn't correct. My two questions are: is my derivation of the gradient correct, and if so, is my implementation in Python correct?  If this is better asked in the Math forum, feel free to move.</p>

<p>Model: <a href="http://i.stack.imgur.com/uyW2a.gif" rel="nofollow"><img src="http://i.stack.imgur.com/uyW2a.gif" alt="enter image description here"></a></p>

<p>Log likelihood: <a href="http://i.stack.imgur.com/SM5wj.gif" rel="nofollow"><img src="http://i.stack.imgur.com/SM5wj.gif" alt="enter image description here"></a></p>

<p>Finally, the gradient:<a href="http://i.stack.imgur.com/TTe3g.gif" rel="nofollow"><img src="http://i.stack.imgur.com/TTe3g.gif" alt="enter image description here"></a></p>

<p>Here, i is each observation, j is an alternative within observation i, c is the chosen alternative in observation i, Xij is the feature vector for choice j in i and B are the corresponding coefficients.  <em>The likelihood formula should have the feature vector multiplied by the coefficient vector.  My mistake</em></p>

<p>My implementation for the likelihood and gradients are as follows: </p>

<p>Likelihood:</p>

<pre><code>def log_likelihood(coefs, observations, config, lasso):
    def func(grp):
        mtrx = grp.as_matrix(config.features)
        dp = np.dot(mtrx, coefs)
        sub = np.log(np.exp(dp).sum())
        inc = (dp * grp['choice']).sum()
        return inc - sub 
    ll = observations.groupby(['observation_id']).apply(func).sum()
    if lasso is not None:
        ll -= (np.abs(coefs).sum() * lasso)
    neg_log = ll * -1
return neg_log
</code></pre>

<p>Gradient:</p>

<pre><code>def gradient(coefs, observations, config, lasso):
    def func(grp):
        mtrx = grp.as_matrix([config.features])
        tmtrx = mtrx.transpose()
        tmp = np.exp(tmtrx * coefs[:, np.newaxis])
        sub = (tmp * tmtrx).sum(1) / tmp.sum(1)
        inc = (mtrx * grp['choice'][:, np.newaxis]).sum(0)
        ret = inc - sub 
        return ret 
    return -1 * observations.groupby(['observation_id']).apply(func).sum()
</code></pre>

<p>Here, coefs is a numpy array containing the coefficients, observations is a dataframe where each row represents an alternative within in an observation and the columns are a choice column indicating 0/1 for the choice within a column and an observation_id column where all alternatives within an observation have the same id, and finally config is a dict containing a member 'features' which is the list of columns in the observations df containing the features. <em>Note I'm testing without using the lasso parameter</em>.  Example below of what the data looks like.</p>

<p>I've verified the likelihood is correct; however, the error for the gradient is very large when using scipy.optimize.check_grad.  I'm also able to solve for B when not passing the gradient to scipy.optimize.minimize.  The gradient evaluates as I'd expect, so at this point I can only think that my derivation is incorrect, but I'm not sure why.</p>

<pre><code>In [27]: df.head(14)
Out[27]: 
          x1        x2        x3  observation_id  choice
0   0.187785  0.435922 -0.475349             211       1
1  -0.935956 -0.405833 -1.753128             211       0
2   0.210424  0.141579  0.415933             211       0
3   0.507025  0.307965 -0.198089             211       0
4   0.080658 -0.125473 -0.592301             211       0
5   0.605302  0.239491  0.287094             293       1
6   0.259580  0.415388 -0.396969             293       0
7  -0.637267 -0.984442 -1.376066             293       0
8   0.241874  0.435922  0.855742             293       0
9   0.831534  0.650425  0.930592             293       0
10 -1.682565  0.435922 -2.517229             293       0
11 -0.149186  0.300299  0.494513             293       0
12 -1.918179 -9.967421 -2.774450             293       0
13 -1.185817  0.295601 -1.974923             293       0
</code></pre>

<p>The derivation was incorrect.  In the exponentiation I was only including the feature and coefficient for the partial derivative of the given coefficient.  Rather, it should have been e to the dot product of all features and coefficients.  </p>

python scipy statistics mathematical-optimization