Density Ratio Trick を使った KL-divergence の計算に関するメモ

##Point
* Density ratio trick による KL-divergence の計算を実装し、シンプルな具体例で確認。

##Sample code
```python
# Density Ratio Trick

def generate_norm_samples(mu, sig, n_samples):
  return np.random.normal(loc = mu, scale = sig, size = (n_samples))

# -----
def kl_norm_analytics(mu_1, sig_1, mu_2, sig_2):
  return np.log(sig_2 /sig_1) + 0.5 * (sig_1**2 + (mu_1 - mu_2)**2 - sig_2**2) / sig_2**2

# -----
def kl_norm_mc(mu_1, sig_1, mu_2, sig_2, n_samples):
  x = generate_norm_samples(mu_1, sig_1, n_samples)
  return np.mean(np.log(pdf_norm(x, mu_1, sig_1 + 1e-8)) \
                        - np.log(pdf_norm(x, mu_2, sig_2) + 1e-8)) 

def pdf_norm(x, mu, sig):
  return 1.0 / np.sqrt(2.0 * np.pi * sig**2) * np.exp(-0.5 * (x - mu)**2 / sig**2)

# -----
def discriminator(x, n_units_1, n_units_2): 
  init_norm = tf.truncated_normal_initializer(stddev = 0.01, dtype = tf.float32)
  with tf.variable_scope('disc'):
    w_1 = tf.get_variable('w_1', shape = [1, n_units_1], initializer = init_norm)
    b_1 = tf.get_variable('b_1', shape = [n_units_1], initializer = init_norm)
    w_2 = tf.get_variable('w_2', shape = [n_units_1, n_units_2], initializer = init_norm)
    b_2 = tf.get_variable('b_2', shape = [n_units_2], initializer = init_norm)
    w_3 = tf.get_variable('w_3', shape = [n_units_2, 1], initializer = init_norm)
    b_3 = tf.get_variable('b_3', shape = [1], initializer = init_norm)
    
    h_1 = tf.nn.tanh(tf.add(tf.matmul(x, w_1), b_1))
    h_2 = tf.nn.tanh(tf.add(tf.matmul(h_1, w_2), b_2))
    h_3 = tf.add(tf.matmul(h_2, w_3), b_3)
  return tf.nn.sigmoid(h_3)

def ce_loss(y_hat, y):
  return - tf.reduce_mean(y * tf.log(tf.clip_by_value(y_hat, 1e-10, 1.0)) \
                          + (1 - y) * tf.log(tf.clip_by_value(1 - y_hat, 1e-10, 1.0)))

def optimize_opt(loss, learning_rate):
  optimizer = tf.train.AdamOptimizer(learning_rate = learning_rate)
  return optimizer.minimize(loss)

def kl_norm_dr(n_units_1, n_units_2, n_iterations, batch_size, show_step, learning_rate, \
          mu_1, sig_1, mu_2, sig_2, n_samples):
  
  tf.reset_default_graph()
  
  x = tf.placeholder(shape = [None, 1], dtype = tf.float32)
  y = tf.placeholder(shape = [None, 1], dtype = tf.float32)
  
  y_hat = discriminator(x, n_units_1, n_units_2)
  loss = ce_loss(y_hat, y)
  opt = optimize_opt(loss, learning_rate)
  
  kl = tf.reduce_mean(tf.log(y_hat + 1e-8) - tf.log(1 - y_hat + 1e-8))
    
  init = tf.global_variables_initializer()
  
  with tf.Session() as sess:
    sess.run(init)
    
    for i in range(n_iterations):
      x_1 = generate_norm_samples(mu_1, sig_1, batch_size // 2)
      x_1 = np.reshape(x_1, (-1, 1)).astype(np.float32)
      y_1 = np.ones(shape = (batch_size // 2, 1), dtype = np.float32)
      x_2 = generate_norm_samples(mu_2, sig_2, batch_size // 2)
      x_2 = np.reshape(x_2, (-1, 1)).astype(np.float32)
      y_2 = np.zeros(shape = (batch_size // 2, 1), dtype = np.float32)

      indices = np.random.choice(batch_size, size = batch_size, replace = False)
      x_batch = np.concatenate([x_1, x_2], axis = 0)[indices]
      y_batch = np.concatenate([y_1, y_2], axis = 0)[indices]
      
      feed_dict = {x: x_batch, y: y_batch}
      
      sess.run(opt, feed_dict = feed_dict)
      
#      if (i + 1) % show_step == 0:
#        temp_loss = sess.run(loss, feed_dict = feed_dict)
#        print ('--------------------')
#        print ('Iteration: ' + str(i + 1) + '  Loss: ' + str(temp_loss))
        
    x_sample = generate_norm_samples(mu_1, sig_1, n_samples)
    x_sample = np.reshape(x_sample, (-1, 1)).astype(np.float32)
    y_hat_sample = sess.run(y_hat, feed_dict = {x: x_sample})
    kl_sample = sess.run(kl, feed_dict = {x: x_sample})

  return kl_sample
```

##Output results
```python
mu_1 = 1.0
sig_1 = 1.0
mu_2 = 2.0
sig_2 = 0.8
n_samples = 500

kl_1 = kl_norm_analytics(mu_1, sig_1, mu_2, sig_2)
print ('Analytics: ' + str(kl_1))

# -----
history_kl = []
for i in range(5):
  x_1 = generate_norm_samples(mu_1, sig_1, n_samples)
  x_2 = generate_norm_samples(mu_2, sig_2, n_samples)

  kl_2 = kl_norm_mc(mu_1, sig_1, mu_2, sig_2, n_samples)
  history_kl.append(kl_2)

kl_2_mean = np.mean(history_kl)
kl_2_std = np.std(history_kl)  

print ('-' * 10)
print ('Monte Carlo mean: ' + str(kl_2_mean))
print ('Monte Carlo std: ' + str(kl_2_std))
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/aad4159a-af36-cdab-e4e4-df58c7eb01e3.png)


```python
n_units_1 = 5
n_units_2 = 5
n_iterations = 60
batch_size = 32
show_step = 20
learning_rate = 0.1

history_kl = []
for i in range(5):
  kl_3 = kl_norm_dr(n_units_1, n_units_2, n_iterations, batch_size, show_step, learning_rate, \
      mu_1, sig_1, mu_2, sig_2, n_samples)
  history_kl.append(kl_3)

kl_3_mean = np.mean(history_kl)
kl_3_std = np.std(history_kl)  

print ('-' * 10)
print ('Density ratio trick mean: ' + str(kl_3_mean))
print ('Density ratio trick std: ' + str(kl_3_std))
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/6964e39d-c284-9cda-13d0-a4340f34f853.png)