计算条件概率分布的累积分布

Question

对于给定的m、p(z|m)，我的条件概率为z，其中选择系数是为了在[0,1.5] 和m 的范围内积分超过z [18:28] 等于一。

def p(z,m):  
    if (m<21.25):
        E = { 'ft':0.55, 'alpha': 2.99, 'z0':0.191, 'km':0.089, 'kt':0.25 }
        S = { 'ft':0.39, 'alpha': 2.15, 'z0':0.121, 'km':0.093, 'kt':-0.175 }
        I={ 'ft':0.06, 'alpha': 1.77, 'z0':0.045, 'km':0.096, 'kt':-0.9196 }
        Evalue=E['ft']*np.exp(-1*E['kt']*(m-18))*z**E['alpha']*np.exp(-1*(z/(E['z0']+E['km']*(m-18)))**E['alpha'])
        Svalue=S['ft']*np.exp(-1*S['kt']*(m-18))*z**S['alpha']*np.exp(-1*(z/(S['z0']+S['km']*(m-18)))**S['alpha'])
        Ivalue=I['ft']*np.exp(-1*I['kt']*(m-18))*z**I['alpha']*np.exp(-1*(z/(I['z0']+I['km']*(m-18)))**I['alpha'])
        value=Evalue+Svalue+Ivalue
    elif(m>=21.25):
        E = { 'ft':0.25, 'alpha': 1.957, 'z0':0.321, 'km':0.196, 'kt':0.565 }
        S = { 'ft':0.61, 'alpha': 1.598, 'z0':0.291, 'km':0.167, 'kt':0.155 }
        I = { 'ft':0.14, 'alpha': 0.964, 'z0':0.170, 'km':0.129, 'kt':0.1759 }
        Evalue=E['ft']*np.exp(-1*E['kt']*(m-18))*z**E['alpha']*np.exp(-1*(z/(E['z0']+E['km']*(m-18)))**E['alpha'])
        Svalue=S['ft']*np.exp(-1*S['kt']*(m-18))*z**S['alpha']*np.exp(-1*(z/(S['z0']+S['km']*(m-18)))**S['alpha'])
        Ivalue=I['ft']*np.exp(-1*I['kt']*(m-18))*z**I['alpha']*np.exp(-1*(z/(I['z0']+I['km']*(m-18)))**I['alpha'])
        value=Evalue+Svalue+Ivalue
    return value

我想从这个分布中抽取一个样本，因此我在z和m平面上做了一个网格点来估计累积分布，m上的累积积分达到1但累积积分超过z 不给我一个在边缘。我不知道为什么它不会收敛到一个？！！

grid_m = np.linspace(18, 28, 1000)
grid_z = np.linspace(0, 1.5, 1000)
dz = np.diff(grid_z[:2])
# get cdf on grid, use cumtrapz 
prob_zgm=np.empty((grid_z.shape[0], grid_m.shape[0]),float)
for i in range(grid_z.shape[0]):
    for j in range(grid_m.shape[0]):
        prob_zgm[i,j]=p(grid_z[i],grid_m[j])

pr = np.column_stack((np.zeros(prob_zgm.shape[0]),prob_zgm))
dm = np.diff(grid_m[:2])
cdf_zgm = integrate.cumtrapz(pr, dx=dm, axis=1)

cdf = integrate.cumtrapz(pr, dx=dz, axis=0)       

哪种假设可能会导致这种不一致或我计算错误？

更新：累积分布cdf_zgm如图所示

剩下的，为了得到概率的倒数，是我用过的方法：

# fix bounds of cdf_zgm
cdf_zgm[:, 0] = 0
cdf_zgm[:, -1] = 1
#Interpolate the data using a linear spline to "grid_q" samples
grid_q = np.linspace(0, 1, 200)
grid_qm = np.empty((len(grid_m), len(grid_q)), float)
for i in range(len(grid_m)):
    grid_qm[i] = interpolate.interp1d(cdf_zgm[i], grid_z)(grid_q)

# build 2d interpolation for z as function of (q,m)
z_interp = interpolate.interp2d(grid_q, grid_m, grid_qm)
#sample magnitude 
ng=20000
r = dist_m.rvs(ng)
rvs_u = np.random.rand(ng)
rvs_z = np.asarray([z_interp(rvs_u[i], r[i]) for i in range(len(rvs_u))]).ravel()

将CDF 的边界固定为one 的方法是否正确？

Answer 1

我不知道该代码有什么问题。但这里有几个不同的想法可以尝试：

(1) 只需对数组元素求和，而不是尝试计算数值积分。这种方式更简单。 （对数组元素求和本质上是计算一个矩形规则近似值，事实证明，它实际上比梯形规则更准确。）

(2) 不要试图一次创建一个完整的二维数组，而是编写一个函数，它只为给定的 m 值创建一个 p(z | m) 的一维切片。然后只需将这些元素相加即可得到累积概率。