How to use the chaospy.E function in chaospy

def test_dist_sub(distribution):
    """Test distribution subtraction."""
    dist1_e = cp.E(distribution() - 3.0)
    dist2_e = cp.E(3.0 - distribution())
    base_e = cp.E(distribution()) - 3.0
    np.testing.assert_allclose(dist1_e, -dist2_e, rtol=1e-05, atol=1e-08)
    np.testing.assert_allclose(dist2_e, -base_e, rtol=1e-05, atol=1e-08)

def test_orth_norms():
    dist = cp.Normal(0, 1)
    orth = cp.orth_ttr(5, dist, normed=True)
    norms = cp.E(orth**2, dist)
    assert np.allclose(norms, 1)

def test_dist_mul(distribution):
    """Test distribution multiplication."""
    dist1_e = cp.E(distribution() * 9.0)
    dist2_e = cp.E(9.0 * distribution())
    base_e = cp.E(distribution()) * 9.0
    np.testing.assert_allclose(dist1_e, dist2_e, rtol=1e-05, atol=1e-08)
    np.testing.assert_allclose(dist2_e, base_e, rtol=1e-05, atol=1e-08)

    dist1_e = cp.E(distribution() * 0.1)
    dist2_e = cp.E(0.1 * distribution())
    base_e = cp.E(distribution()) * 0.1
    np.testing.assert_allclose(dist1_e, dist2_e, rtol=1e-05, atol=1e-08)
    np.testing.assert_allclose(dist2_e, base_e, rtol=1e-05, atol=1e-08)

def test_operator_E():
    dist = cp.Normal(0, 1)
    res = np.array([1, 0, 1, 0, 3])
    x = cp.variable()
    poly = x**np.arange(5)
    assert np.allclose(cp.E(poly, dist), res)

plt.fill_between(
        coord, expected-deviation, expected+deviation,
        color="k", alpha=0.3
    )
    plt.plot(coord, expected, "k--", lw=3)
    plt.xlabel(r"\verb;coord;")
    plt.ylabel(r"function evaluations \verb;foo;")
    plt.title("Results using point collocation method")
    plt.savefig("results_collocation.png")
    plt.clf()


    absissas, weights = cp.generate_quadrature(8, distribution, "C")
    evals = [foo(coord, val) for val in absissas.T]
    foo_approx = cp.fit_quadrature(polynomial_expansion, absissas, weights, evals)
    expected = cp.E(foo_approx, distribution)
    deviation = cp.Std(foo_approx, distribution)

    plt.fill_between(
        coord, expected-deviation, expected+deviation,
        color="k", alpha=0.3
    )
    plt.plot(coord, expected, "k--", lw=3)
    plt.xlabel(r"\verb;coord;")
    plt.ylabel(r"function evaluations \verb;foo;")
    plt.title("Results using psuedo-spectral projection method")
    plt.savefig("results_spectral.png")
    plt.clf()

P_nk = cp.outer(P, P)
E_ank = cp.E(q0*P_nk, dist)
E_ik = cp.E(q1*P, dist)
sE_ank = cp.sum(E_ank, 0)

# Right hand side of the ODE
def f(c_k, x):
    return -cp.sum(c_k*E_ank, -1)/norm

solver = odespy.RK4(f)
c_0 = E_ik/norm
solver.set_initial_condition(c_0)
c_n, x = solver.solve(x)
U_hat = cp.sum(P*c_n, -1)

mean = cp.E(U_hat, dist)
var = cp.Var(U_hat, dist)

# actually start the postprocessing now:

nodes, weights = cp.generate_quadrature(order, distribution, rule='G',sparse=sparse)
expansion,norms = cp.orth_ttr(3, distribution,retall=True)
approx_denit = cp.fit_quadrature(expansion, nodes, weights, np.mean(data[:,1,:], axis=1))
approx_oxy = cp.fit_quadrature(expansion, nodes, weights, np.mean(data[:,0,:], axis=1))

annual_oxy = cp.fit_quadrature(expansion,nodes,weights,data[:,0,:])
annual_denit = cp.fit_quadrature(expansion,nodes,weights,data[:,1,:])

s_denit = cp.descriptives.sensitivity.Sens_m(annual_denit, distribution)
s_oxy = cp.descriptives.sensitivity.Sens_m(annual_oxy, distribution)

df_oxy = cp.Std(annual_oxy,distribution)
df_denit = cp.Std(annual_denit,distribution)
f0_oxy = cp.E(annual_oxy,distribution)
f0_denit = cp.E(annual_denit,distribution)

#  = <-a*sum(c*P), P[k]>
# d/dx c[k]* = -sum(c*)
# d/dx c = -E( outer(a*P,P) ) / E( P*P )
#
# u(0) = I
#  = 
# c[k](0)  = 
# c(0) = E( I*P ) / E( P*P )
order = 5
P, norm = cp.orth_ttr(order, dist, retall=True, normed=True)

# support structures
q0, q1 = cp.variable(2)
P_nk = cp.outer(P, P)
E_ank = cp.E(q0*P_nk, dist)
E_ik = cp.E(q1*P, dist)
sE_ank = cp.sum(E_ank, 0)

# Right hand side of the ODE
def f(c_k, x):
    return -cp.sum(c_k*E_ank, -1)/norm

solver = odespy.RK4(f)
c_0 = E_ik/norm
solver.set_initial_condition(c_0)
c_n, x = solver.solve(x)
U_hat = cp.sum(P*c_n, -1)

mean = cp.E(U_hat, dist)
var = cp.Var(U_hat, dist)