Mean Prediction Interval Width (MPIW) and Kullback-Leibler Divergence (KL-Divergence)

[Estol111]

Hi, author. When I calculate the evaluation metrics, I am very confused about the determination of the upper bound U and the lower bound L. Would it be convenient for you to provide the codes for MPIW and KL-Divergence? Thank you!

[ZhuangDingyi]

Upper and lower bounds are calculated through sampling. For example, after you read the data:

from scipy import stats

stzinb_res = np.load('../output/ny_full_60min_rand0_ZISTNB_nonorm.npz')
n_stzinb = stzinb_res['n']
p_stzinb = stzinb_res['p']
pi_stzinb = stzinb_res['pi']
maxval_stzinb = stzinb_res['max_value']
tar_stzinb = stzinb_res['target']
pred_stzinb = stzinb_res['mean_pred']

pi_idx1 = pi_stzinb[s:e,od,0]<=0.01
pi_idx2 = (pi_stzinb[s:e,od,0]>0.01) & (pi_stzinb[s:e,od,0]<0.99) 

L_stzinb = np.zeros_like(pi_stzinb[s:e,od,0])
U_stzinb = np.zeros_like(pi_stzinb[s:e,od,0])

U_stzinb[pi_idx1] = stats.nbinom.ppf(0.01-pi_stzinb[s:e,od,0][pi_idx1], n_stzinb[s:e,od,0][pi_idx1], p_stzinb[s:e,od,0][pi_idx1])
U_stzinb[pi_idx2] = 0
L_stzinb[pi_idx1] = stats.nbinom.ppf(0.99-pi_stzinb[s:e,od,0][pi_idx1], n_stzinb[s:e,od,0][pi_idx1], p_stzinb[s:e,od,0][pi_idx1])
L_stzinb[pi_idx2] = stats.nbinom.ppf(0.99-pi_stzinb[s:e,od,0][pi_idx2], n_stzinb[s:e,od,0][pi_idx2], p_stzinb[s:e,od,0][pi_idx2])

L_stnb = np.zeros_like(n_stnb[s:e,od,0])
U_stnb = np.zeros_like(n_stnb[s:e,od,0])

U_stnb = stats.nbinom.ppf(0.01, n_stnb[s:e,od,0], p_stnb[s:e,od,0])
L_stnb = stats.nbinom.ppf(0.99, n_stnb[s:e,od,0], p_stnb[s:e,od,0])

I would not recommend using KL-Divergence as I realize it is not suitable for the evaluation of point distributions in our case, but MPIW is fine. I define the KL-Divergence in the whole dataset level, which is deviated from the original definition of it. But here is the code block I used for the paper:

def KL_DIV_divide(truth,pred): return np.sum( pred*np.log( (pred+1e-5)/(truth+1e-5) ) )/np.prod(truth.shape)

[Estol111]

Thank you very much！

[Estol111]

When I tried to apply the code you provided to my task, I found that when calculating the upper and lower bounds, the output was 0. I suspect this may have something to do with the three parameters of the ZINB model I used that were predicted by the deep learning model. I'm wondering if I need to do an anti-normalization process to get the right result? If you have ever been in a similar situation or have any suggestions, I would appreciate your guidance.

[ZhuangDingyi]

Not always, if you train correctly. This is the code of plotting Figure 4(b):

fig,ax = universal_fig_row_col(2,1,figsize=(2.3,3),fontsize=6,sharex=False,sharey=False) # change figsize
ax[0].scatter(Demand,MPIW_stzinb,label='STZINB',s=0.2,color='C1')
ax[0].scatter(Demand,MPIW_stnb,label='STNB',s=0.2,color='C4')
ax[0].scatter(Demand,MPIW_stgau,label='STG',s=0.2,color='C0')
ax[0].set_yscale('symlog')
ax[0].scatter(Demand,MPIW_sttruncnorm,label='STTN',s=0.2,color='C3')
ax[0].legend(markerscale=10.,loc=4,frameon=True)
ax[0].set_yticklabels(['0','1','10'])
ax[0].set_xlabel('Average trip demand per O-D (trips/15min)')
ax[0].set_ylabel('MPIW')
# ax[0].set_title('SLD_15min')

s = 3000 + 19*4 + 4#410;1200 # Apr 25
e = 3300 #510;2600
od = 1325
lw = 1.5
ax[1].plot(tar_stnb[s:e,od,0],label='True',linewidth=lw,color='C0')
ax[1].plot(pred_stzinb[s:e,od,0].astype(np.int),label='STZINB',linewidth=lw,color='C1')
ax[1].plot(pred_median_stnb[s:e,od,0].astype(np.int),label='STNB',linewidth=lw,color='C4')

pi_idx1 = pi_stzinb[s:e,od,0]<=0.01
pi_idx2 = (pi_stzinb[s:e,od,0]>0.01) & (pi_stzinb[s:e,od,0]<0.99) 
L_stzinb = np.zeros_like(pi_stzinb[s:e,od,0])
U_stzinb = np.zeros_like(pi_stzinb[s:e,od,0])

U_stzinb[pi_idx1] = stats.nbinom.ppf(0.01-pi_stzinb[s:e,od,0][pi_idx1], n_stzinb[s:e,od,0][pi_idx1], p_stzinb[s:e,od,0][pi_idx1])
U_stzinb[pi_idx2] = 0
L_stzinb[pi_idx1] = stats.nbinom.ppf(0.99-pi_stzinb[s:e,od,0][pi_idx1], n_stzinb[s:e,od,0][pi_idx1], p_stzinb[s:e,od,0][pi_idx1])
L_stzinb[pi_idx2] = stats.nbinom.ppf(0.99-pi_stzinb[s:e,od,0][pi_idx2], n_stzinb[s:e,od,0][pi_idx2], p_stzinb[s:e,od,0][pi_idx2])

L_stnb = np.zeros_like(n_stnb[s:e,od,0])
U_stnb = np.zeros_like(n_stnb[s:e,od,0])

U_stnb = stats.nbinom.ppf(0.01, n_stnb[s:e,od,0], p_stnb[s:e,od,0])
L_stnb = stats.nbinom.ppf(0.99, n_stnb[s:e,od,0], p_stnb[s:e,od,0])

plt.fill_between(np.arange(0,e-s),U_stnb.astype(np.int),L_stnb.astype(np.int),color='C4',alpha=0.3)
plt.fill_between(np.arange(0,e-s),U_stzinb.astype(np.int),L_stzinb.astype(np.int),color='C1',alpha=0.3)

ax[1].legend(loc=1,frameon=True)
# ax[1].set_xlabel('t/15 min (will change into dates)')
ax[1].set_xticks(np.arange(0,e-s,4*8))
ax[1].set_xticklabels(['0:00\nApr 25','8:00','16:00','0:00\nApr 26','8:00','16:00','0:00\nApr 27'])
ax[1].set_ylabel('Trip demand')
plt.tight_layout()
plt.savefig('fig/uq_15min.pdf')

I have no normalization process as you need to fit your data directly to calculate the likelihood. Some points might be fitted with both 0 upper and lower bounds, but not always the case. If your data are too sparse, it is likely that your lower bounds are always 0, but the upper bound should be some value larger than 0.

[Estol111]

Thank you very much! I find the problem!