Make CoxPHFitter linear instead of quadratic

This commit optimizes the algorithm which calculates Efron gradients and
hessian during the Newton-Raphson fitting.

The change basically comes down to noting the fact that:

    R(t(i-1)) = [t(i-1)] + R(t(i))

where (the last risk set) R(t(n)) = [t(n)]

The math and original implementation translates to something like:

    # Nested sums as in original equations: O(n^2) running time
    for i in range(n):
        phi_sum = 0
        for j in riskset(i):
            phi_sum += phi(j)

Whereas taking advantage of the recursive relationship is akin to
something like:

    # Use recursive relationship: O(n) running time
    phi_sum = 0
    for i in reversed(range(n)):
        phi_sum += phi(i)

The key is to iterate over the data backwards. For this to work, the
data must be sorted on time (O(nlog(n))). I have compared the old and
new implementations on a real data set (2695,12) and they give the same
results, ignoring normal floating point differences due to order of
summation.

Performance-wise, for data with shape (2695,12), the fitting speed goes
from 6min 27sec to 1.3sec on my machine (9 iterations of N-R).

I am still a bit unfamiliar with Pandas and therefor had some indexing troubles.
That's why the arrays are forced to numpy arrays in _newton_. This could
perhaps be improved.
lol so embarrassing, thanks for fixing that.