Geek asks and answers

How to solve a one-parameter equation using Python (scipy / numpy?)

I hope you have some useful advice so that I come up with the following task:

I wrote some simple python snippets to build probability density functions. In my particular case, let them represent conditionally conditional probabilities for some parameter x.

So, I am wondering if there is a smart approach (i.e. a module) in Python (possibly through a function or method of NumPy or SciPy) to solve a simple equation for a parameter x. For instance,

pdf(x, mu=10, sigma=3**0.5) / pdf(x, mu=20, sigma=2**0.5) = 1
# get x

Right now, I can only use brute force when I use something like x = np.arange(0, 50, 0.000001)and save the value of x in a vector that gives the closest value to 1 when calculating the ratiopdf1/pdf2.

Below is the code that I wrote to calculate pdf and built the relation:

def pdf(x, mu=0, sigma=1):
    """Calculates the normal distribution probability density 
        function (PDF).  

    """
    term1 = 1.0 / ( math.sqrt(2*np.pi) * sigma )
    term2 = np.exp( -0.5 * ( (x-mu)/sigma )**2 )
    return term1 * term2


x = np.arange(0, 100, 0.05)

pdf1 = pdf(x, mu=10, sigma=3**0.5)
pdf2 = pdf(x, mu=20, sigma=2**0.5)

# ...
# ratio = pdf1 / pdf2
# plt.plot(x, ratio)

Thank!

0
source share
2 answers

In general, it looks like you need scalar root lookup functions: http://docs.scipy.org/doc/scipy/reference/optimize.html

But, as others have noted, it seems that there is an analytical solution.

+1
source

Since you have a good closed-form equation, you can solve it with SymPy.

I hooked up the values ​​for muand sigmaand injected this into Sympy Gamma :

 solve(1.0 / ( sqrt(2*pi) *(3**0.5) ) * exp( -0.5 * ( (x-10)/(3**0.5) )**2 ) /  (1.0 / ( sqrt(2*pi) *(2**0.5) ) * exp( -0.5 * ( (x-20)/(2**0.5) )**2 ))-1,x)

Result: 15.4554936768195

+1
source

More articles:


All Articles