我正在尝试对 1/H(z)dz 进行数值积分，其中函数 H(z) 由下式给出
def Hz_th(theta,z):
aa,H0,omegam0=theta
return np.sqrt(H0**2*(1 + z)**3*omegam0 + 3*H0**2*(1 + z)**3*aa*omegam0 + (1 + z)**3*np.sqrt((4*H0**4*(1 - omegam0))/(1 + z)**6 + (-(H0**2*omegam0) - 3*H0**2*aa*omegam0)**2))/np.sqrt(2)
其中 z 是一个数组
z = np.array([0.50349, 0.4952, 0.6782, 0.8672, 0.7992, 0.37129, 0.35568, 0.28391, 0.46691, 0.94791, 0.69391, 0.89791, 0.78991, 0.59091, 0.2102, 0.8592, 0.5592, 0.9972, 0.5834, 0.7972, 0.9142, 0.94921, 0.7342, 0.7662, 0.9842, 0.7012, 0.83716, 0.37016, 0.35816, 0.52216, 0.70116, 0.74526, 0.59116, 0.62039, 1.2024, 0.6104, 0.35819, 0.45138, 0.55239, 0.75039, 0.82238, 0.95039, 0.34039, 0.8104, 0.61191, 0.54891, 1.02991, 0.42291, 0.70191, 0.68591, 0.93291, 0.74991, 0.46891, 0.95991, 0.8412, 0.5632, 0.8292, 0.57921, 0.8652, 0.8592, 0.58621, 0.4892, 0.68921, 0.32416, 0.18516, 0.68116, 0.92116, 0.42816, 0.69916, 0.57516, 0.51116, 0.64116, 0.35116, 0.34916, 0.35486, 0.87116, 0.92516, 0.4194, 0.89038, 0.48038, 0.58038, 0.80538, 0.66438, 0.74538, 0.7364, 0.37039, 0.85039, 0.21938, 0.64739, 0.60038, 0.76039, 0.47039, 0.49791, 0.50791, 0.58791, 0.69992, 0.76191, 0.78891, 0.40591, 0.80892, 0.58392, 0.72091, 0.83791, 0.80691, 0.50016, 0.53316, 0.82116, 0.12516, 0.44316, 0.68339, 0.72639, 0.40439, 0.69038, 0.57639, 0.72039, 0.55091, 0.8489, 0.43591, 0.30191, 0.58119, 0.34621, 0.3312, 0.44939, 0.2916, 0.46109, 0.63291, 0.26891, 0.92591, 0.62591, 0.60891, 0.57891, 0.7202, 0.7672, 0.3682, 0.8492, 0.77721, 0.58421, 0.5892, 0.7692, 0.5142, 0.63821, 0.43521, 0.58421, 0.66259, 0.62116, 0.41616, 0.45117, 0.73316, 0.74116, 0.34916, 0.51416, 0.64339, 0.4704, 0.26338, 0.7424, 0.91039, 0.3374, 0.98339, 0.81338, 0.8174, 0.96038, 0.75639, 0.87891, 0.80991, 0.62791, 0.93491, 0.69791, 0.76791, 0.86491, 0.51491, 0.47091, 0.92891, 0.6312, 0.6162, 0.5652, 0.55819, 0.7362, 0.2622, 0.7622, 0.8592, 0.5592, 0.48016, 0.92116, 0.47516, 0.89216, 0.63116, 0.73516, 0.41816, 0.53516, 0.73416, 0.60916, 0.70116, 0.98216, 0.75816, 0.64339, 0.51539, 0.80538, 0.80039, 0.71538, 0.76639, 0.90139, 0.71839, 0.57938, 0.96038, 0.64839, 0.93639, 0.46138, 0.67039, 0.24639, 0.63892, 0.53591, 0.77391, 0.37091, 0.69992, 0.37391, 0.85391, 0.60392, 0.53491, 0.40091, 0.37191, 0.84091, 0.53216, 0.93116, 0.55316, 0.84116, 0.72739, 0.44238, 0.2824, 0.5194, 1.06039, 0.26838, 0.2504, 0.69892, 0.73091, 0.99891, 0.40991, 0.12694, 0.25609, 0.08848, 0.11797, 0.18322, 0.1407, 0.15186, 0.09401, 0.28662, 0.26349, 0.1431, 0.16035, 0.12303, 0.12608, 0.17315, 0.16431, 0.1599, 0.24448, 0.24835, 0.27544, 0.33051, 0.33103, 0.16059, 0.21887, 0.11991, 0.17839, 0.27923, 0.18581, 0.06641, 0.258, 0.2798, 0.26576, 0.185, 0.05715, 0.3264, 0.10886, 0.20323, 0.37986, 0.2116, 0.21353, 0.14409, 0.25249, 0.25672, 0.21885, 0.07883, 0.08282, 0.17028, 0.2367, 0.30569, 0.08254, 0.18891, 0.08141, 0.25655, 0.12622, 0.10264, 0.12455, 0.37046, 0.20215, 0.1799, 0.1618, 0.22669, 0.29683, 0.25029, 0.14882, 0.14729, 0.21513, 0.14478, 0.13243, 0.21975, 0.38492, 0.174, 0.17789, 0.29936, 0.11845, 0.17653, 0.24552, 0.19611, 0.24852, 0.30632, 0.24637, 0.13514, 0.21054, 0.2363, 0.17851, 0.22068, 0.31494, 0.18443, 0.26129, 0.18123, 0.38048, 0.18481, 0.20847, 0.17706, 0.31582, 0.24973, 0.20368, 0.12336, 0.20318, 0.15339, 0.33439, 0.1958, 0.24973, 0.15509, 0.15596, 0.15837, 0.30938, 0.24857, 0.18507, 0.18376, 0.18079, 0.17696, 0.17768, 0.0718, 0.13811, 0.1696, 0.24555, 0.14387, 0.25191, 0.03743, 0.27749, 0.28746, 0.19867, 0.15999, 0.21075, 0.23775, 0.15517, 0.09391, 0.15271, 0.25722, 0.21226, 0.2454, 0.13696, 0.32569, 0.10638, 0.40127, 0.18931, 0.19726, 0.12837, 0.12262, 0.12715, 0.21323, 0.21099, 0.26856, 0.20453, 0.2701, 0.24185, 0.29114, 0.35346, 0.28543, 0.21955, 0.17558, 0.23542, 0.15325, 0.21005, 0.29135, 0.31376, 0.13687, 0.29567, 0.30394, 0.20557, 0.1219, 0.17428, 0.24883, 0.24558, 0.13833, 0.24315, 0.24267, 0.31468, 0.27573, 0.30041, 0.21728, 0.20335, 0.19025, 0.29021, 0.22786, 0.31065, 0.13729, 0.08784, 0.08504, 0.19009, 0.09227, 0.12903, 0.26166, 0.26162, 0.19422, 0.14703, 0.21179, 0.17958, 0.19029, 0.3388, 0.11741, 0.2889, 0.26405, 0.24961, 0.22967, 0.08543, 0.27656, 0.36193, 0.30021, 0.11635, 0.15492, 0.25037, 0.25174, 0.12928, 0.30929, 0.29888, 0.27091, 0.29353, 0.18979, 0.12376, 0.31288, 0.30915, 0.20267, 0.21163, 0.1789, 0.30354, 0.21457, 0.32045, 0.21835, 0.07523, 0.2576, 0.18468, 0.11628, 0.06481, 0.10288, 0.04437, 0.16595, 0.16817, 0.17063, 0.12042, 0.13868, 0.2338, 0.24643, 0.22286, 0.17915, 0.19718, 0.22916, 0.26451, 0.10337, 0.21338, 0.15653, 0.2434, 0.23652, 0.20852, 0.25104, 0.22534, 0.28742, 0.32984, 0.24603, 0.11966, 0.32931, 0.21082, 0.20308, 0.28436, 0.18037, 0.2607, 0.14792, 0.32793, 0.24206, 0.22334, 0.1828, 0.15304, 0.13818, 0.20735, 0.20313, 0.19801, 0.22256, 0.2734, 0.18651, 0.15935, 0.21935, 0.17373, 0.24436, 0.13353, 0.28332, 0.17383, 0.19654, 0.36952, 0.20548, 0.17357, 0.12688, 0.28434, 0.19392, 0.10018, 0.37425, 0.11771, 0.19176, 0.23684, 0.15899, 0.07877, 0.13864, 0.23381, 0.15689, 0.08806, 0.13045, 0.18203, 0.25625, 0.21674, 0.14533, 0.13639, 0.20515, 0.36976, 0.1294, 0.26687, 0.27822, 0.17483, 0.11364, 0.11242, 0.15317, 0.19671, 0.18914, 0.17625, 0.21757, 0.18563, 0.16356, 0.25721, 0.12164, 0.16455, 0.19862, 0.26821, 0.40128, 0.14581, 0.05573, 0.22452, 0.28655, 0.18388, 0.10351, 0.33136, 0.21865, 0.26357, 0.16284, 0.10683, 0.16477, 0.1075, 0.12899, 0.05948, 0.04093, 0.01705, 0.01531, 0.01472, 0.02673, 0.01732, 0.04437, 0.01567, 0.02477, 0.038, 0.0299, 0.01012, 0.0258, 0.01038, 0.03791, 0.03507, 0.01268, 0.02525, 0.01908, 0.03489, 0.0295, 0.03416, 0.04067, 0.07515, 0.05581, 0.02529, 0.03501, 0.02846, 0.0328, 0.03378, 0.04597, 0.04341, 0.05684, 0.02024, 0.03191, 0.0415, 0.02639, 0.02651, 0.02689, 0.02349, 0.05067, 0.04093, 0.06877, 0.01754, 0.02557, 0.02189, 0.03155, 0.06845, 0.02249, 0.04949, 0.03161, 0.03434, 0.04486, 0.01531, 0.03048, 0.03837, 0.03682, 0.0221, 0.01402, 0.0589, 0.05895, 0.06888, 0.0204, 0.02894, 0.02384, 0.01451, 0.03292, 0.06455, 0.01663, 0.05359, 0.03163, 0.02751, 0.02559, 0.02401, 0.03713, 0.03408, 0.02454, 0.01697, 0.01629, 0.03184, 0.037, 0.02932, 0.0174, 0.03578, 0.02845, 0.05035, 0.04269, 0.03212, 0.02752, 0.03049, 0.01724, 0.07063, 0.05344, 0.02878, 0.03518, 0.03348, 0.04651, 0.06119, 0.04714, 0.0212, 0.04028, 0.03863, 0.03314, 0.02414, 0.02837, 0.01043, 0.02038, 0.0471, 0.03984, 0.0122, 0.01271, 0.01763, 0.02171, 0.01572, 0.01972, 0.01678, 0.01226, 0.05322, 0.03055, 0.03001, 0.01506, 0.0153, 0.04713, 0.03411, 0.0152, 0.01504, 0.02763, 0.03359, 0.01447, 0.02547, 0.01082, 0.01741, 0.02031, 0.04819, 0.02223, 0.02811, 0.05458, 0.05721, 0.02091, 0.01454, 0.0431, 0.03273, 0.03399, 0.01382, 0.05092, 0.01752, 0.02482, 0.0288, 0.01568, 0.02197, 0.03426, 0.01698, 0.02752, 0.03438, 0.01318, 0.01726, 0.01209, 0.05418, 0.0197, 0.01662, 0.02398, 0.03144, 0.01291, 0.03827, 0.01532, 0.02643, 0.0202, 0.0373, 0.02348, 0.01635, 0.02275, 0.2305, 0.3812, 0.2534, 0.1526, 0.4365, 0.1023, 0.3712, 0.4661, 0.1959, 0.3112, 0.0322, 0.1098, 0.2255, 0.2604, 0.2253, 0.4474, 0.4399, 0.3374, 0.5283, 0.5761, 0.3876, 0.36, 0.347, 0.4391, 0.3492, 0.4305, 0.5088, 0.2888, 0.1727, 0.318, 0.2752, 0.103, 0.2049, 0.2647, 0.2348, 0.2145, 0.4276, 0.4019, 0.3218, 0.4512, 0.1612, 0.3512, 0.2492, 0.0822, 0.1056, 0.3327, 0.1022, 0.2928, 0.1326, 0.2707, 0.1004, 0.3536, 0.4235, 0.3001, 0.1476, 0.2288, 0.1388, 0.1436, 0.0919, 0.1568, 0.2388, 0.1496, 0.5078, 0.2528, 0.2842, 0.2368, 0.3412, 0.2641, 0.349, 0.3212, 0.1352, 0.3065, 0.2335, 0.3006, 0.3676, 0.1692, 0.07976, 0.1806, 0.1806, 0.2307, 0.04582, 0.2212, 0.2006, 0.0755, 0.02627, 0.1912, 0.1777, 0.1852, 0.3491, 0.2712, 0.2921, 0.3704, 0.0963, 0.1808, 0.2748, 0.2454, 0.3618, 0.2318, 0.1466, 0.1592, 0.4198, 0.3192, 0.5196, 0.1506, 0.3026, 0.2106, 0.2396, 0.3316, 0.2006, 0.3606, 0.2426, 0.2212, 0.1206, 0.3412, 0.2504, 0.3804, 0.3106, 0.1994, 0.24, 0.5027, 0.3701, 0.27, 0.2631, 0.268, 0.308, 0.25, 0.2804, 0.3148, 0.3288, 0.3297, 0.2588, 0.2246, 0.2102, 0.3878, 0.3068, 0.3652, 0.2488, 0.2988, 0.1238, 0.3312, 0.3826, 0.2406, 0.1712, 0.2088, 0.2896, 0.2457, 0.1361, 0.1505, 0.1187, 0.2448, 0.3503, 0.6192, 0.2196, 0.0994, 0.4098, 0.5025, 0.3815, 0.2797, 0.2307, 0.4269, 0.5501, 0.3253, 0.3509, 0.3283, 0.4504, 0.06441, 0.06559, 0.3185, 0.1453, 0.1396, 0.4294, 0.2166, 0.3666, 0.3806, 0.4206, 0.1996, 0.5106, 0.4267, 0.3156, 0.5779, 0.4812, 0.3922, 0.1735, 0.3302, 0.2996, 0.2296, 0.2927, 0.2085, 0.2749, 0.3032, 0.359, 0.36, 0.349, 0.6304, 0.3207, 0.51, 0.3669, 0.1384, 0.1792, 0.1717, 0.2708, 0.3473, 0.4796, 0.3192, 0.5192, 0.3802, 0.3712, 0.1808, 0.2236, 0.1412, 0.3055, 0.4045, 0.2856, 0.3012, 0.2637, 0.5006, 0.3051, 0.3857, 0.1985, 0.0906, 0.4774, 0.3312, 0.4712, 0.3306, 0.2607, 0.07141, 0.1812, 0.3379, 0.1221, 0.2501, 0.1604, 0.2988, 0.3892, 0.2928, 0.3288, 0.3056, 0.2858, 0.4292, 0.1792, 0.1896, 0.2488, 0.4096, 0.2416, 0.3376, 0.0706, 0.3456, 0.3306, 0.05495, 0.4231, 0.1072, 0.2854, 0.3254, 0.3347, 0.2907, 0.4404, 0.2411, 0.4607, 0.1194, 0.545, 0.4114, 0.2688, 0.04976, 0.4118, 0.07369, 0.3188, 0.4198, 0.4208, 0.4392, 0.4804, 0.1998, 0.3018, 0.2902, 0.3092, 0.2956, 0.3386, 0.3782, 0.0896, 0.4507, 0.4604, 1.206, 1.33, 1.54, 1.55, 1.7, 1.8, 2.26, 1.914, 1.3, 1.34, 01.02, 0.735, 1.12, 1.23, 1.23, 0.854, 1.37, 0.9752, 0.97, 0.74, 1.39, 1.305, 0.935, 1.014, 1.315, 1.092])
此外，我需要使用积分对 θ 中包含的自由参数进行 MCMC 采样。我尝试了以下代码来对数组执行数值积分
def MCMCf(theta,z):
aa,H0,omegam0=theta
return np.vectorize(lambda x: 1.0/Hz_th(theta,z))
def MCMCfint(theta,z):
aa,H0,omegam0=theta
return np.vectorize(integrate.quad(MCMCf(theta,z), 0, z))
但这对我不起作用。也许有人可以提供建议？我需要保留函数形式，因为我无法提供自由参数的值，它们将受到
emcee当我尝试用随机真值作为自由参数绘制结果时，有一个完整的回溯
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/tmp/ipykernel_17393/1708689059.py in <module>
7 return np.vectorize(integrate.quad(MCMCf(aa,H0,omegam0,z), 0, z))
8 z=z_values
----> 9 plt.plot(z,MCMCfint(aa_true,H0_true,omegam0_true,z))
/tmp/ipykernel_17393/1708689059.py in MCMCfint(aa, H0, omegam0, z)
5
6 def MCMCfint(aa,H0,omegam0,z):
----> 7 return np.vectorize(integrate.quad(MCMCf(aa,H0,omegam0,z), 0, z))
8 z=z_values
9 plt.plot(z,MCMCfint(aa_true,H0_true,omegam0_true,z))
~/.local/lib/python3.8/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
346
347 # check the limits of integration: \int_a^b, expect a < b
--> 348 flip, a, b = b < a, min(a, b), max(a, b)
349
350 if weight is None:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
import scipy.special as spec
from scipy.integrate
import quad
import numpy as np
array = np.array([0, 1, 2, 3, 4, 5]) # array[i] >= 0
for all i in [0, len(array) - 1].For example, array = np.array([1, 2, 3, 4, 5])
constant = 3 # >= 0. For example, constant = 3
a = array
b = constant * a # multiply numpy array by a scalar(it is possible with numpy
arrays(: ) c = b * a / constant
e3 = spec.expn(3, c)
print('Output:', e3) #
return a ndarray with the value of the integral e3 en each point
specified in c(numpy array)
# Output: [1.09691967e-01 2.76136095e-03 1.04786279e-05 5.96855753e-09 4.97790975e-13]
# For evaluate the integral e3 from 0 to a[len(a) - 1]:
result1 = e3[len(a) - 1] - e3[0] print('Output:', result1) # Output: -0.499999999950223
# Other thing is the integral of e3 from 0 to a[len(a) - 1].You can use
quad(
return the value and an estimted error) # See https:
result2 = quad(lambda x: spec.expn(3, x), 0, c[len(a) - 1])
print('Output:', result2) # Output: (0.3333333333328521, 4.440892098500626e-16)