Solving Power Law Distribution in Python

I have data that closely resembles a power law distribution. Using Python, I want to approximate the data by solving two equations in the form:

y is the y axis data. In Python it would be data[i]. x would be i + 1. It follows that we get two equations with two unknown variables at the first data index and at a "random" 2nd one somewhere else in the data:

The problem comes down to solving just

due to mathematical simplification. I don't know how to solve an equation like this using libraries like numpy.linalg.solve. How do I find the value of a using Python?

Solution

Alright, I got it.

import math

def get_power_law_variables(data):
    c = data[0]
    middle_index = len(data) / 2
    division = float(data[middle_index]) / c
    logarithm_base = middle_index + 1
    a = math.log(division, logarithm_base)
    return c, a

# Example usage
data = range(50, 150)
c, a = get_power_law_variables(data)
print c, a