Как рассчитать корреляцию между двумя распределениями в байесовской статистике? ⇐ Python

[Anonymous]

Я пытаюсь вычислить корреляцию между двумя распределениями (назовем их ростом и весом). Я знаю, как вычислить простую корреляцию Пирсона, используя частотный подход: вычислить ковариацию, вычислить ковариационную матрицу, а затем разделить ее на sd(x)×sd(y)
Стандартную формулу ковариации можно рассчитать как: Формула ковариации 1
или через коэффициент корреляции 𝜌 как:
Формула ковариации 2

• Первый вопрос, который я Вопрос в том, уместно ли рассчитывать среднее значение, как в частотном методе. Было бы правильно получить два отдельных апостериорных показателя (для веса и роста), а затем рассчитать ковариацию и корреляцию, используя первое уравнение?
• Правильно ли это? можно ли рассматривать ρ как апостериорный? Должен ли я сначала получить апостериорные данные о росте и весе, а затем уже рассчитывать ρ?

Мой код >
Я пытаюсь использовать PyMC, но в моей модели я могу получить оценку только для одного дистрибутива. Априорный и апостериорный ρ различаются, но я не понимаю, почему.

Код: Выделить всё

high_values = df['high'].values #no NaN
weight_values = df["weight"].values

with pm.Model() as model:
rho_transformed = pm.Beta("rho_transformed", alpha=2, beta=2)
rho = pm.Deterministic("rho", rho_transformed * 2 - 1)  # Map back to [-1, 1]

# Define mean
mu_high = pm.NegativeBinomial("mu_1", mu=mu_1, alpha = 4) #t or poisson ...
mu_weight = pm.NegativeBinomial('mu_2', mu=mu_2, alpha=shape) #negBin
# Difine std
sigma_high = pm.HalfNormal("sigma_1", sigma=1)
sigma_weight = pm.HalfNormal("sigma_2", sigma=1)

cov_matrix = pm.math.stack([
[sigma_high**2, rho * sigma_high * sigma_weight],
[rho * sigma_high * sigma_weight, sigma_weight**2]
]) #is it correct in bayes?

# Observed data
observed_data = np.column_stack([high_values, weight_values]).astype("float64")

mu_value = pt.as_tensor_variable([mu_high, mu_weight]) #thats a very strange thing
observed = pm.MvNormal("observed", mu=mu_value, cov=cov_matrix, observed=observed_data) #multivariate normal distribution

# Prior and Posterior Predictive Checks
prior_samples = pm.sample_prior_predictive(draws=1000)
posterior_samples = pm.sample_posterior_predictive(trace)

# Run the sampler
trace = pm.sample(1000, return_inferencedata=True) #Markov Chain Monte Carlo (MCMC) sampling to draw samples from the posterior distribution of each parameter.

Вы можете вызвать Prior_samples.prior['rho'] и трассировку.posterior['rho'], чтобы проверить, отличаются ли они:)

Подробнее здесь: https://stackoverflow.com/questions/791 ... statistics