ValueError when using scipy.integrate.quad

Question

I am receiving a ValueError when integrating using scipy.integrate.quad. Here is my simplified code:

import numpy as np
import scipy.integrate as integrate
p = np.arange(0,1,1/1000)
h = lambda p: p**2/2+p*(1-p)
Kl = lambda p: h(p) + 0.02
K = Kl(p)
R = 0.5*h(p) + 0.5*h(1)
Vl = lambda p: np.minimum.reduce([p, K, R])
integrate.quad(Vl, 0, 1)[0]

The last line gives the exception:

ValueError: setting an array element with a sequence.

Can someone please propose how to do this integration properly? Thanks

Answer 1

I think this code confuses functions on floats with arrays. The function argument to quad needs to accept a float and return a float. In

Vl = lambda p: np.minimum.reduce([p, K, R])

p is a float (and is not the p array-- it's usually a bad idea to have a module global with the same name as an argument) and K and R are arrays, which isn't what we want (that's what's throwing the error.)

I think you're just looking for

def h(p):
    return p**2/2+p*(1-p)

def K(p):
    return h(p) + 0.02

def R(p):
    return 0.5*h(p) + 0.5*h(1)

def Vl(p):
    return min(p, K(p), R(p))

which gives me

In [177]: integrate.quad(Vl, 0, 1.0)
Out[177]: (0.34689543041336846, 4.8736376714649885e-09)