0

For the purpose of my problem, I have to use gamma functions from mpmath module rather than scipy. However, I am using matplotlib.pyplot to plot my functions when the variable is taken from a numpy array and the y-values are supposed to be produced in the same format. But in doing so, I am receiving an error message of TypeError complaining of the conversion from array to mpf. 

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad, dblquad
import mpmath as mp


low, up  = 5.630e5, 1.167e12
alpha, threshold = 1.05   , 2.15e10     
beta = 274

def g(x, low, up, beta):
    return mp.gamma(-2/3) * (mp.gammainc(-2/3, beta*(x/low)**3) - mp.gammainc(-2/3, beta*(x/up)**3))

def Integrand1(x, low, up, threshold, alpha, beta):
    return pow(x/threshold, alpha) * g(x, low, up, beta)

def Integrand2(x, low, up, threshold, alpha, beta):
    return g(x, low, up, beta)


def PDF(x, low, up, threshold, alpha, beta):
    A = quad(Integrand1, low, threshold, args=(low, up, threshold, alpha, beta))
    B = quad(Integrand2, threshold, up, args=(low, up, threshold, alpha, beta))
    C=(A[0]+B[0])**(-1)
    y = np.piecewise(x, 
                    [x < threshold], [lambda x: C * pow(x/threshold, alpha) * g(x, low, up, beta), 
                                   lambda x: C *                           g(x, low, up, beta)
                                  ]
                    )
    return y


x_array = np.array(np.logspace(8.2, 11.8, 10))
y_array = PDF(x_array, low, up, threshold, alpha, beta)
plt.plot(x_array, y_array, color='green', linestyle='-')
plt.gca().autoscale(False)
plt.vlines([threshold], 
        plt.gca().get_ylim()[0], plt.gca().get_ylim()[1], 
        linestyles='dashed', color='k', label='')

plt.xscale("log", nonposx='clip')
plt.yscale("log", nonposy='clip')
plt.show()

Traceback (most recent call last): File "test.py", line 35, in y_array = PDF(x_array, low, up, threshold, alpha, beta) File "test.py", line 28, in PDF lambda x: C * g(x, low, up, beta) File "/home/username/anaconda3/lib/python3.6/site-packages/numpy/lib/function_base.py", line 1344, in piecewise y[condlist[k]] = item(vals, *args, **kw) File "test.py", line 27, in [x < threshold], [lambda x: C * pow(x/threshold, alpha) * g(x, low, up, beta), File "test.py", line 13, in g return mp.gamma(-2/3) * (mp.gammainc(-2/3, beta*(x/low)**3) - mp.gammainc(-2/3, beta*(x/up)**3)) File "/home/username/anaconda3/lib/python3.6/site-packages/mpmath/functions/expintegrals.py", line 141, in gammainc a = ctx.convert(a) File "/home/username/anaconda3/lib/python3.6/site-packages/mpmath/ctx_mp_python.py", line 662, in convert return ctx._convert_fallback(x, strings) File "/home/username/anaconda3/lib/python3.6/site-packages/mpmath/ctx_mp.py", line 614, in _convert_fallback raise TypeError("cannot create mpf from " + repr(x)) TypeError: cannot create mpf from array([ 6.11259157e+09, 9.68780477e+10,
1.53541358e+12, 2.43346654e+13, 3.85678455e+14, 6.11259157e+15])

I need to know how to mpf an array and I found thes pages (How to mpf an array) and (create an mpf) but I don't know how to apply them to my plotting routine.

Rebel
  • 472
  • 8
  • 25
  • After remove the function *g* and replacing it simply with unity, the code works. So obviously, something has to be done for the function g so that it is compatible with the other functions. – Rebel Nov 21 '17 at 05:41

1 Answers1

1

This isn't about plotting. The error occurs within the function PDF, specifically within np.piecewise. As the documentation says, the lambdas you pass to piecewise must be able to take an array, because this is what piecewise will pass to them.

And your function g isn't able to take an array, because the mpmath functions it uses are meant for a single float. Solution: np.vectorize wrapper:

gv = np.vectorize(g)

and then use gv instead of g in piecewise.

Then np.piecewise will work, and you'll be able to proceed to other errors, like mismatched names PDF_values and y_array.

  • @thanks. After correcting that vecorization, I am receiving this: "TypeError: unsupported operand type(s) for *: 'float' and 'vectorize'" in the piecewise line. All I did replaced g with gv and the statement you mentioned right after the definition of the g function. – Rebel Nov 21 '17 at 07:33
  • Works on my computer: `y = np.piecewise(x, [x < threshold], [lambda x: C * pow(x/threshold, alpha) * gv(x, low, up, beta), lambda x: C * gv(x, low, up, beta) ] )` –  Nov 21 '17 at 07:37
  • Explicit variables of gv were missing in the parentheses. I thought it could memorize it from the defined statement. – Rebel Nov 21 '17 at 07:44