Генерация и оценка 200 нормалей и их уменьшение - PullRequest
Новая
       93

Генерация и оценка 200 нормалей и их уменьшение

0 голосов
/ 29 апреля 2020

Я пытаюсь оценить нормальную плотность, используя квадратичное c приближение в тензорном потоке (код 4.14 из Статистического переосмысления McElreath).

Код, который у меня есть на данный момент: 

import pandas as pd
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp
from  tensorflow_probability import distributions as tfd

_BASE_URL = "https://raw.githubusercontent.com/rmcelreath/rethinking/Experimental/data"

 HOWELL_DATASET_PATH = f"{_BASE_URL}/Howell1.csv"

df = pd.read_csv(HOWELL_DATASET_PATH, sep=';')
df = df[df['age'] >= 18]

mu = tf.linspace(start=140.0, stop=160.0, num=200)
sigma= tf.linspace(start=4.0, stop=9.0, num=200)

tf.reduce_sum(tfd.Normal(loc=mu, scale=sigma).log_prob(df.height))

Сбой из-за того, что df имеет форму (352,), в то время как я создаю (200,) точек для моего нормального распределения в

Однако 

tf.reduce_sum(tfd.Normal(loc=mu, scale=sigma).log_prob(2))

и 

tf.reduce_sum(tfd.Normal(loc=mu[0], scale=sigma[0]).log_prob(df.height))

оба работают.

Мне нужно создать (200, 352) тензор - один нормальный для каждого mu, sigma в моей сетке, а затем оценить его с моими данными выборки - df. У меня вопрос: как мне это сделать?

Ответы [ 2 ]

1 голос
/ 01 мая 2020

Я думаю, что совместное распространение TFP - хороший способ express это: 

mu = tf.linspace(start=140.0, stop=160.0, num=200)
sigma = tf.linspace(start=7.0, stop=9.0, num=200)

def mk_joint(nobs):
  return tfd.JointDistributionNamed(dict(
      mu=tfd.Normal(178, 20),
      sigma=tfd.Uniform(0, 50),
      height=lambda mu, sigma: tfd.Sample(tfd.Normal(loc=mu, scale=sigma), nobs)
  ))
joint = mk_joint(len(df))
joint.sample()
print(f'joint event shape: {joint.event_shape}')
lp = joint.log_prob(dict(mu=mu[:,tf.newaxis], sigma=sigma, height=df.height))
import matplotlib.pyplot as plt
plt.imshow(lp)
plt.xlabel('sigma')
plt.xticks(np.arange(len(sigma))[::10], sigma[::10].numpy().round(2), rotation=90)
plt.ylabel('mu')
plt.yticks(np.arange(len(mu))[::10], mu[::10].numpy().round(2))
plt.show()

=> joint event shape: {'sigma': TensorShape([]), 'mu': TensorShape([]), 'height': TensorShape([352])}

enter image description here

0 голосов
/ 29 апреля 2020

Итак, я выяснил, что один из способов сделать это - создать сетку (200, 200, 352), а затем изменить ее форму, а остальные вычисления будут выполнены напрямую. 

import pandas as pd
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp
from  tensorflow_probability import distributions as tfd

_BASE_URL = "https://raw.githubusercontent.com/rmcelreath/rethinking/Experimental/data"

 HOWELL_DATASET_PATH = f"{_BASE_URL}/Howell1.csv"

df = pd.read_csv(HOWELL_DATASET_PATH, sep=';')
df = df[df['age'] >= 18]


mu = tf.linspace(start=140.0, stop=160.0, num=200)
sigma = tf.linspace(start=7.0, stop=9.0, num=200)

means, variances, _  = tf.meshgrid(mu, sigma,  np.zeros((352,)).astype(np.float32))
means = tf.reshape(means, [40000, 352])
variances = tf.reshape(variances, [40000, 352])

normal = tfd.Normal(loc=means, scale=variances)

log_lik = tf.reduce_sum(normal.log_prob(df.height), axis=1)

logprob_mu = tfd.Normal(178.0, 20.0).log_prob(means)
logprob_sigma = tfd.Uniform(low=0.0, high=50.0).log_prob(variances)

log_joint_prod = log_lik + logprob_mu[:, 0] + logprob_sigma[:, 0]
joint_prob_tf = tf.exp(log_joint_prod - tf.reduce_max(log_joint_prod))
...