Antes de entrarmos nos detalhes técnicos: este submit é, obviamente, dedicado a McElreath, que escreveu um dos livros mais intrigantes sobre a modelagem bayesiana (ou devemos apenas dizer – científica?) Da qual estamos cientes. Se você não leu Repensar estatísticae estão interessados em modelar, você pode definitivamente querer conferir. Neste submit, não vamos tentar rever a história: nosso foco claro será, em vez disso, uma demonstração de como fazer o MCMC com TFProbability.
Concretamente, este submit tem duas partes. O primeiro é uma visão geral rápida de como usar TFD_JONT_SESTEAL_DISTRIBUTION Para construir um modelo e depois provar com ele usando o Hamiltoniano Monte Carlo. Esta parte pode ser consultada para pesquisa rápida de código ou um modelo frugal de todo o processo. A segunda parte passa por um modelo de vários níveis com mais detalhes, mostrando como extrair, pós-processo e visualizar amostragem, bem como saídas de diagnóstico.
ReedFrogs
Os dados vêm com o rethinking pacote.
'information.body': 48 obs. of 5 variables:
$ density : int 10 10 10 10 10 10 10 10 10 10 ...
$ pred : Issue w/ 2 ranges "no","pred": 1 1 1 1 1 1 1 1 2 2 ...
$ measurement : Issue w/ 2 ranges "huge","small": 1 1 1 1 2 2 2 2 1 1 ...
$ surv : int 9 10 7 10 9 9 10 9 4 9 ...
$ propsurv: num 0.9 1 0.7 1 0.9 0.9 1 0.9 0.4 0.9 ...A tarefa é modelar a contagem de sobreviventes entre os girinos, onde os girinos são mantidos em tanques de tamanhos diferentes (equivalentemente, diferentes números de habitantes). Cada linha no conjunto de dados descreve um tanque, com sua contagem inicial de habitantes (density) e número de sobreviventes (surv). Na parte de visão geral técnica, construímos um modelo simples não -frio que descreve todos os tanques isolados. Então, no passo a passo detalhado, veremos como construir um interceptações variadas modelo que permite o compartilhamento de informações entre os tanques.
Construindo modelos com tfd_joint_distribution_sequential
tfd_joint_distribution_sequential representa um modelo como uma lista de distribuições condicionais. É mais fácil ver em um exemplo actual, por isso entraremos direto, criando um modelo não frio dos dados de girinos.
É assim que a especificação do modelo seria Stan:
mannequin{
vector(48) p;
a ~ regular( 0 , 1.5 );
for ( i in 1:48 ) {
p(i) = a(tank(i));
p(i) = inv_logit(p(i));
}
S ~ binomial( N , p );
}E aqui está tfd_joint_distribution_sequential:
library(tensorflow)
# be sure you have not less than model 0.7 of TensorFlow Likelihood
# as of this writing, it's required of set up the grasp department:
# install_tensorflow(model = "nightly")
library(tfprobability)
n_tadpole_tanks <- nrow(d)
n_surviving <- d$surv
n_start <- d$density
m1 <- tfd_joint_distribution_sequential(
record(
# regular prior of per-tank logits
tfd_multivariate_normal_diag(
loc = rep(0, n_tadpole_tanks),
scale_identity_multiplier = 1.5),
# binomial distribution of survival counts
operate(l)
tfd_independent(
tfd_binomial(total_count = n_start, logits = l),
reinterpreted_batch_ndims = 1
)
)
)O modelo consiste em duas distribuições: meios anteriores e variações para os 48 tanques de girotos são especificados por tfd_multivariate_normal_diag; então tfd_binomial gera contagens de sobrevivência para cada tanque. Observe como a primeira distribuição é incondicional, enquanto a segunda depende do primeiro. Observe também como o segundo deve ser envolvido tfd_independent Para evitar transmissão errada. (Este é um aspecto de tfd_joint_distribution_sequential O uso que merece ser documentado mais sistematicamente, o que certamente vai acontecer. Apenas pense que essa funcionalidade foi adicionada ao TFP grasp Apenas três semanas atrás!)
Como um aparte, a especificação do modelo aqui acaba mais curta do que em Stan como tfd_binomial Opcionalmente, leva logits como parâmetros.
Como em cada distribuição de TFP, você pode fazer uma verificação rápida de funcionalidade, amostragem do modelo:
# pattern a batch of two values
# we get samples for each distribution within the mannequin
s <- m1 %>% tfd_sample(2)((1))
Tensor("MultivariateNormalDiag/pattern/affine_linear_operator/ahead/add:0",
form=(2, 48), dtype=float32)
((2))
Tensor("IndependentJointDistributionSequential/pattern/Beta/pattern/Reshape:0",
form=(2, 48), dtype=float32)e probabilidades de log de computação:
# we must always get solely the general log likelihood of the mannequin
m1 %>% tfd_log_prob(s)t((1))
Tensor("MultivariateNormalDiag/pattern/affine_linear_operator/ahead/add:0",
form=(2, 48), dtype=float32)
((2))
Tensor("IndependentJointDistributionSequential/pattern/Beta/pattern/Reshape:0",
form=(2, 48), dtype=float32)Agora, vamos ver como podemos provar deste modelo usando o Hamiltonian Monte Carlo.
Correndo Hamiltonian Monte Carlo em TFP
Definimos um kernel Hamiltoniano Monte Carlo com adaptação dinâmica de tamanho de etapa com base em uma probabilidade de aceitação desejada.
# variety of steps to run burnin
n_burnin <- 500
# optimization goal is the chance of the logits given the information
logprob <- operate(l)
m1 %>% tfd_log_prob(record(l, n_surviving))
hmc <- mcmc_hamiltonian_monte_carlo(
target_log_prob_fn = logprob,
num_leapfrog_steps = 3,
step_size = 0.1,
) %>%
mcmc_simple_step_size_adaptation(
target_accept_prob = 0.8,
num_adaptation_steps = n_burnin
)Em seguida, executamos o amostrador, passando em um estado inicial. Se queremos correr (n ) correntes, esse estado tem que ter de comprimento (n )para todos os parâmetros do modelo (aqui temos apenas um).
A função de amostragem, mcmc_sample_chainpode opcionalmente ser passado um trace_fn Isso informa ao TFP quais tipos de meta informações a salvar. Aqui, economizamos índices de aceitação e tamanhos de etapas.
# variety of steps after burnin
n_steps <- 500
# variety of chains
n_chain <- 4
# get beginning values for the parameters
# their form implicitly determines the variety of chains we'll run
# see current_state parameter handed to mcmc_sample_chain beneath
c(initial_logits, .) %<-% (m1 %>% tfd_sample(n_chain))
# inform TFP to maintain monitor of acceptance ratio and step measurement
trace_fn <- operate(state, pkr) {
record(pkr$inner_results$is_accepted,
pkr$inner_results$accepted_results$step_size)
}
res <- hmc %>% mcmc_sample_chain(
num_results = n_steps,
num_burnin_steps = n_burnin,
current_state = initial_logits,
trace_fn = trace_fn
)Quando a amostragem é concluída, podemos acessar as amostras como res$all_states:
mcmc_trace <- res$all_states
mcmc_traceTensor("mcmc_sample_chain/trace_scan/TensorArrayStack/TensorArrayGatherV3:0",
form=(500, 4, 48), dtype=float32)Esta é a forma das amostras para los 48 logits por tanques: 500 amostras vezes 4 cadeias vezes 48 parâmetros.
A partir dessas amostras, podemos calcular o tamanho eficaz da amostra e (rhat ) (Alias mcmc_potential_scale_reduction):
# Tensor("Imply:0", form=(48,), dtype=float32)
ess <- mcmc_effective_sample_size(mcmc_trace) %>% tf$reduce_mean(axis = 0L)
# Tensor("potential_scale_reduction/potential_scale_reduction_single_state/sub_1:0", form=(48,), dtype=float32)
rhat <- mcmc_potential_scale_reduction(mcmc_trace)Enquanto que as informações de diagnóstico estão disponíveis em res$hint:
# Tensor("mcmc_sample_chain/trace_scan/TensorArrayStack_1/TensorArrayGatherV3:0",
# form=(500, 4), dtype=bool)
is_accepted <- res$hint((1))
# Tensor("mcmc_sample_chain/trace_scan/TensorArrayStack_2/TensorArrayGatherV3:0",
# form=(500,), dtype=float32)
step_size <- res$hint((2)) Após esse esboço rápido, vamos para o tópico prometido no título: modelagem de vários níveis ou pool parcial. Desta vez, também examinaremos mais de perto os resultados da amostra e saídas de diagnóstico.
Girinos de vários níveis
O modelo de vários níveis-ou Modelo de interceptações variadasneste caso: vamos chegar a inclinações variadas Em um submit posterior – adiciona um Hiperprior para o modelo. Em vez de decidir sobre uma média e variação do regular, as logits são desenhadas, deixamos o modelo aprender meios e variações para tanques individuais. Esses meios por tanques, apesar de serem anteriores aos logits binomiais, assumem-se que são normalmente distribuídos e são regularizados por um antes da média e um exponencial antes da variação.
Para o Stan-Eavvy, aqui está a formulação STAN deste modelo.
mannequin{
vector(48) p;
sigma ~ exponential( 1 );
a_bar ~ regular( 0 , 1.5 );
a ~ regular( a_bar , sigma );
for ( i in 1:48 ) {
p(i) = a(tank(i));
p(i) = inv_logit(p(i));
}
S ~ binomial( N , p );
}E aqui está com TFP:
m2 <- tfd_joint_distribution_sequential(
record(
# a_bar, the prior for the imply of the conventional distribution of per-tank logits
tfd_normal(loc = 0, scale = 1.5),
# sigma, the prior for the variance of the conventional distribution of per-tank logits
tfd_exponential(price = 1),
# regular distribution of per-tank logits
# parameters sigma and a_bar consult with the outputs of the above two distributions
operate(sigma, a_bar)
tfd_sample_distribution(
tfd_normal(loc = a_bar, scale = sigma),
sample_shape = record(n_tadpole_tanks)
),
# binomial distribution of survival counts
# parameter l refers back to the output of the conventional distribution instantly above
operate(l)
tfd_independent(
tfd_binomial(total_count = n_start, logits = l),
reinterpreted_batch_ndims = 1
)
)
)Tecnicamente, dependências em tfd_joint_distribution_sequential são definidos por proximidade espacial na lista: no aprendizado anterior para os logits
operate(sigma, a_bar)
tfd_sample_distribution(
tfd_normal(loc = a_bar, scale = sigma),
sample_shape = record(n_tadpole_tanks)
)sigma refere -se à distribuição imediatamente acima e a_bar para o acima disso.
Analogamente, na distribuição das contagens de sobrevivência
operate(l)
tfd_independent(
tfd_binomial(total_count = n_start, logits = l),
reinterpreted_batch_ndims = 1
)l refere -se à distribuição imediatamente anterior à sua própria definição.
Novamente, vamos provar deste modelo para ver se as formas estão corretas.
s <- m2 %>% tfd_sample(2)
s Eles são.
((1))
Tensor("Regular/sample_1/Reshape:0", form=(2,), dtype=float32)
((2))
Tensor("Exponential/sample_1/Reshape:0", form=(2,), dtype=float32)
((3))
Tensor("SampleJointDistributionSequential/sample_1/Regular/pattern/Reshape:0",
form=(2, 48), dtype=float32)
((4))
Tensor("IndependentJointDistributionSequential/sample_1/Beta/pattern/Reshape:0",
form=(2, 48), dtype=float32)E para garantir que obtemos um no geral log_prob por lote:
Tensor("JointDistributionSequential/log_prob/add_3:0", form=(2,), dtype=float32)Treinamento Este modelo funciona como antes, exceto que agora o estado inicial compreende três parâmetros, a_barAssim, Sigma e l:
c(initial_a, initial_s, initial_logits, .) %<-% (m2 %>% tfd_sample(n_chain))Aqui está a rotina de amostragem:
# the joint log likelihood now's primarily based on three parameters
logprob <- operate(a, s, l)
m2 %>% tfd_log_prob(record(a, s, l, n_surviving))
hmc <- mcmc_hamiltonian_monte_carlo(
target_log_prob_fn = logprob,
num_leapfrog_steps = 3,
# one step measurement for every parameter
step_size = record(0.1, 0.1, 0.1),
) %>%
mcmc_simple_step_size_adaptation(target_accept_prob = 0.8,
num_adaptation_steps = n_burnin)
run_mcmc <- operate(kernel) {
kernel %>% mcmc_sample_chain(
num_results = n_steps,
num_burnin_steps = n_burnin,
current_state = record(initial_a, tf$ones_like(initial_s), initial_logits),
trace_fn = trace_fn
)
}
res <- hmc %>% run_mcmc()
mcmc_trace <- res$all_statesDesta vez, mcmc_trace é uma lista de três: temos
((1))
Tensor("mcmc_sample_chain/trace_scan/TensorArrayStack/TensorArrayGatherV3:0",
form=(500, 4), dtype=float32)
((2))
Tensor("mcmc_sample_chain/trace_scan/TensorArrayStack_1/TensorArrayGatherV3:0",
form=(500, 4), dtype=float32)
((3))
Tensor("mcmc_sample_chain/trace_scan/TensorArrayStack_2/TensorArrayGatherV3:0",
form=(500, 4, 48), dtype=float32)Agora, vamos criar nós gráficos para os resultados e informações em que estamos interessados.
# as above, that is the uncooked end result
mcmc_trace_ <- res$all_states
# we carry out some reshaping operations instantly in tensorflow
all_samples_ <-
tf$concat(
record(
mcmc_trace_((1)) %>% tf$expand_dims(axis = -1L),
mcmc_trace_((2)) %>% tf$expand_dims(axis = -1L),
mcmc_trace_((3))
),
axis = -1L
) %>%
tf$reshape(record(2000L, 50L))
# diagnostics, additionally as above
is_accepted_ <- res$hint((1))
step_size_ <- res$hint((2))
# efficient pattern measurement
# once more we use tensorflow to get conveniently formed outputs
ess_ <- mcmc_effective_sample_size(mcmc_trace)
ess_ <- tf$concat(
record(
ess_((1)) %>% tf$expand_dims(axis = -1L),
ess_((2)) %>% tf$expand_dims(axis = -1L),
ess_((3))
),
axis = -1L
)
# rhat, conveniently post-processed
rhat_ <- mcmc_potential_scale_reduction(mcmc_trace)
rhat_ <- tf$concat(
record(
rhat_((1)) %>% tf$expand_dims(axis = -1L),
rhat_((2)) %>% tf$expand_dims(axis = -1L),
rhat_((3))
),
axis = -1L
) E estamos prontos para realmente executar as correntes.
# to this point, no sampling has been completed!
# the precise sampling occurs after we create a Session
# and run the above-defined nodes
sess <- tf$Session()
eval <- operate(...) sess$run(record(...))
c(mcmc_trace, all_samples, is_accepted, step_size, ess, rhat) %<-%
eval(mcmc_trace_, all_samples_, is_accepted_, step_size_, ess_, rhat_)Desta vez, vamos realmente inspecionar esses resultados.
Girinos de vários níveis: resultados
Primeiro, como as correntes se comportam?
Lotes de rastreamento
Extrair as amostras para a_bar e sigmabem como um dos anteriores instruídos para os logits:
Aqui está um enredo de traço para a_bar:
prep_tibble <- operate(samples) {
as_tibble(samples, .name_repair = ~ c("chain_1", "chain_2", "chain_3", "chain_4")) %>%
add_column(pattern = 1:500) %>%
collect(key = "chain", worth = "worth", -pattern)
}
plot_trace <- operate(samples, param_name) {
prep_tibble(samples) %>%
ggplot(aes(x = pattern, y = worth, colour = chain)) +
geom_line() +
ggtitle(param_name)
}
plot_trace(a_bar, "a_bar")
E aqui para sigma e a_1:
Que tal as distribuições posteriores dos parâmetros, em primeiro lugar, as interceptações variadas a_1 … a_48?
Distribuições posteriores
plot_posterior <- operate(samples) {
prep_tibble(samples) %>%
ggplot(aes(x = worth, colour = chain)) +
geom_density() +
theme_classic() +
theme(legend.place = "none",
axis.title = element_blank(),
axis.textual content = element_blank(),
axis.ticks = element_blank())
}
plot_posteriors <- operate(sample_array, num_params) {
plots <- purrr::map(1:num_params, ~ plot_posterior(sample_array( , , .x) %>% as.matrix()))
do.name(grid.organize, plots)
}
plot_posteriors(mcmc_trace((3)), dim(mcmc_trace((3)))(3))
Agora, vamos ver as médias posteriores correspondentes e os mais altos intervalos de densidade posterior. (O código abaixo inclui os hiperPRIORs em abstract Como queremos exibir um completo summary-como a saída em breve.)
Meios posteriores e hpdis
all_samples <- all_samples %>%
as_tibble(.name_repair = ~ c("a_bar", "sigma", paste0("a_", 1:48)))
means <- all_samples %>%
summarise_all(record (~ imply)) %>%
collect(key = "key", worth = "imply")
sds <- all_samples %>%
summarise_all(record (~ sd)) %>%
collect(key = "key", worth = "sd")
hpdis <-
all_samples %>%
summarise_all(record(~ record(hdi(.) %>% t() %>% as_tibble()))) %>%
unnest()
hpdis_lower <- hpdis %>% choose(-incorporates("higher")) %>%
rename(lower0 = decrease) %>%
collect(key = "key", worth = "decrease") %>%
organize(as.integer(str_sub(key, 6))) %>%
mutate(key = c("a_bar", "sigma", paste0("a_", 1:48)))
hpdis_upper <- hpdis %>% choose(-incorporates("decrease")) %>%
rename(upper0 = higher) %>%
collect(key = "key", worth = "higher") %>%
organize(as.integer(str_sub(key, 6))) %>%
mutate(key = c("a_bar", "sigma", paste0("a_", 1:48)))
abstract <- means %>%
inner_join(sds, by = "key") %>%
inner_join(hpdis_lower, by = "key") %>%
inner_join(hpdis_upper, by = "key")
abstract %>%
filter(!key %in% c("a_bar", "sigma")) %>%
mutate(key_fct = issue(key, ranges = distinctive(key))) %>%
ggplot(aes(x = key_fct, y = imply, ymin = decrease, ymax = higher)) +
geom_pointrange() +
coord_flip() +
xlab("") + ylab("submit. imply and HPDI") +
theme_minimal() 
Agora para um equivalente a summary. Já calculamos meios, desvios padrão e o intervalo HPDI. Vamos adicionar n_effo número efetivo de amostras e rhata estatística de Gelman-Rubin.
Resumo abrangente (também conhecido como “summary”)
is_accepted <- is_accepted %>% as.integer() %>% imply()
step_size <- purrr::map(step_size, imply)
ess <- apply(ess, 2, imply)
summary_with_diag <- abstract %>% add_column(ess = ess, rhat = rhat)
summary_with_diag# A tibble: 50 x 7
key imply sd decrease higher ess rhat
1 a_bar 1.35 0.266 0.792 1.87 405. 1.00
2 sigma 1.64 0.218 1.23 2.05 83.6 1.00
3 a_1 2.14 0.887 0.451 3.92 33.5 1.04
4 a_2 3.16 1.13 1.09 5.48 23.7 1.03
5 a_3 1.01 0.698 -0.333 2.31 65.2 1.02
6 a_4 3.02 1.04 1.06 5.05 31.1 1.03
7 a_5 2.11 0.843 0.625 3.88 49.0 1.05
8 a_6 2.06 0.904 0.496 3.87 39.8 1.03
9 a_7 3.20 1.27 1.11 6.12 14.2 1.02
10 a_8 2.21 0.894 0.623 4.18 44.7 1.04
# ... with 40 extra rows Para as interceptações variadas, o tamanho eficaz de amostras é bastante baixo, indicando que podemos querer investigar possíveis motivos.
Vamos também exibir probabilidades de sobrevivência posterior, analogamente à Figura 13.2 no livro.
Probabilidades de sobrevivência posterior
sim_tanks <- rnorm(8000, a_bar, sigma)
tibble(x = sim_tanks) %>% ggplot(aes(x = x)) + geom_density() + xlab("distribution of per-tank logits")
# our ordinary sigmoid by one other title (undo the logit)
logistic <- operate(x) 1/(1 + exp(-x))
probs <- map_dbl(sim_tanks, logistic)
tibble(x = probs) %>% ggplot(aes(x = x)) + geom_density() + xlab("likelihood of survival")
Finalmente, queremos ter certeza de ver o comportamento de encolhimento exibido na Figura 13.1 no livro.
Encolhimento
abstract %>%
filter(!key %in% c("a_bar", "sigma")) %>%
choose(key, imply) %>%
mutate(est_survival = logistic(imply)) %>%
add_column(act_survival = d$propsurv) %>%
choose(-imply) %>%
collect(key = "kind", worth = "worth", -key) %>%
ggplot(aes(x = key, y = worth, colour = kind)) +
geom_point() +
geom_hline(yintercept = imply(d$propsurv), measurement = 0.5, colour = "cyan" ) +
xlab("") +
ylab("") +
theme_minimal() +
theme(axis.textual content.x = element_blank())
Vemos resultados semelhantes em espírito às de McElreath: as estimativas são encolhidas para a média (a linha cor de ciano). Além disso, o encolhimento parece ser mais ativo em tanques menores, que são os de menor número à esquerda da trama.
Panorama
Neste submit, vimos como construir um interceptações variadas modelo com tfprobabilitybem como como extrair resultados de amostragem e diagnósticos relevantes. Em um próximo submit, seguiremos para inclinações variadas. Com probabilidade não negligenciável, nosso exemplo se baseará em um dos MC Elreath’s novamente … obrigado pela leitura!