1

In python, I would like to find the roots of equations of the form:

-x*log(x) + (1-x)*log(n) - (1-x)*log(1 - x) - k = 0

where n and k are parameters that will be specified.

An additional constraint on the roots is that x >= (1-x)/n. So just for what it's worth, I'll be filtering out roots that don't satisfy that.

My first attempt was to use scipy.optimize.fsolve (note that I'm just setting k and n to be 0 and 1 respectively):

def f(x):                                                                       
    return -x*log(x) + (1-x)*log(1) - (1-x)*log(1-x)                            

fsolve(f, 1)

Using math.log, I got value-errors because I was supplying bad input to log. Using numpy.log gave me some divide by zeros and invalid values in multiply.

I adjusted f as so, just to see what it would do:

def f(x):                                                                       
    if x <= 0:                                                                  
        return 1000                                                             
    if x >= 1:                                                                  
        return 2000                                                             
    return -x*log(x) + (1-x)*log(1) - (1-x)*log(1-x)

Now I get

/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:221: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)

Using python, how can I solve for x for various n and k parameters in the original equation?

oadams
  • 3,019
  • 6
  • 30
  • 53
  • 1
    You have positivity constraints because of the log. I'd change variables for y = log (x / (1 - x)) since x here must be between 0 and 1. (look for logit function on google) – Alexandre C. Dec 03 '13 at 00:02

1 Answers1

1

fsolve also allows guesses to be inserted for where to start. My suggestion would be to plot the equation and have the user type a initial guess either with the mouse or via text to use as an initial guess. You may also want to change the out of bounds values:

if x <= 0:
    return 1000 + abs(x)
if x >= 1:
    return 2000 + abs(x)

This way the function has a slope outside of the region of interest that will guide the solver back into the interesting region.

Troy Rockwood
  • 5,875
  • 3
  • 15
  • 20
  • Hi Troy. There won't be any users interacting with the software. Yeah, I'm having trouble determining what sort of values to use as starting points. 0.5 yields 0.48309238, 0.75 yields 0.99999997. To my understanding these values of n and k should result in x being 0 or 1. – oadams Dec 03 '13 at 00:38
  • hmm, I'm unsure where the 0.48... number comes from, it's not even the max of the function (which is around 0.559). I'm unsure what method fsolve uses in order to find the roots. It may be easier to write your own implementation or use someone else's: http://stackoverflow.com/questions/13054758/python-finding-multiple-roots-of-nonlinear-equation – Troy Rockwood Dec 03 '13 at 00:55