# Bounding JD in terms of TV

a guest
Dec 26th, 2020
78
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import sys, re, math, os, json, random
2. import numpy as np
3.
4. def tvd(x, y):
5.     return 0.5 * sum(abs(np.array(x) - np.array(y)))
6.
7. def jpd(x, y):
8.     ret = 0.0
9.     for i in range(len(x)):
10.         if x[i] == 0 or y[i] == 0:
11.             continue
12.         ret0 = 0.0
13.         for j in range(len(x)):
14.             ret0 += max(x[j] / x[i], y[j] / y[i])
15.         ret += 1.0 / ret0
16.     return 1.0 - ret
17.
18.
19. def simplex_sample(n):
20.     x = np.random.exponential(scale=10.0, size=n)
21.     return x / sum(x)
22.
23.
24. record1 = -100
25. record2 = 100
26. for _ in range(1000000):
27.     x = simplex_sample(3)
28.     y = simplex_sample(3)
29.     tv = tvd(x, y)
30.     jp = jpd(x, y)
31.     ratio1 = jp / (tv / (1 + tv))
32.     ratio2 = jp / tv
33.     if ratio1 > record1:
34.         print('upper bound: %0.10f' % ratio1)
35.         record1 = ratio1
36.     if ratio2 < record2:
37.         print('lower bound: %0.10f' % ratio2)
38.         record2 = ratio2