Em um Postagem anteriormostramos como usar TFProbability – a interface R à probabilidade de tensorflow – para construir um Multinívelou agrupamento parcial Modelo de sobrevivência de girinos em tanques de tamanho diferente (e, portanto, diferentes em número habitante).
Um completamente reunido o modelo teria resultado em uma estimativa international da contagem de sobrevivência, independentemente do tanque, enquanto um não -frio O modelo teria aprendido a prever a contagem de sobrevivência para cada tanque separadamente. A abordagem anterior não leva em consideração diferentes circunstâncias; Este último não utiliza informações comuns. (Além disso, ele claramente não tem uso preditivo, a menos que desejemos fazer previsões para as mesmas entidades que usamos para treinar o modelo.)
Por outro lado, a parcialmente agrupado O modelo permite fazer previsões para as entidades familiares e novas: basta usar o anterior apropriado.
Assumindo que nós são De fato, interessado nas mesmas entidades – por que queremos aplicar um pool parcial? Pelas mesmas razões, tanto esforço no aprendizado de máquina entra para a concepção de mecanismos de regularização. Não queremos demais para ajustar muito para medições reais, elas estão relacionadas à mesma entidade ou a uma classe de entidades. Se eu quiser prever minha frequência cardíaca enquanto acordo na manhã seguinte, com base em uma única medição que estou tomando agora (digamos que é noite e estou digitando freneticamente uma postagem no weblog), é melhor levar em consideração alguns fatos sobre o comportamento da freqüência cardíaca em geral (em vez de apenas projetar no futuro o valor exato medido agora).
No exemplo do girino, isso significa que esperamos que a generalização funcione melhor para tanques com muitos habitantes, em comparação com ambientes mais solitários. Para os últimos, é melhor dar uma olhada nas taxas de sobrevivência de outros tanques, para complementar as informações escassas e idiossincráticas disponíveis. Ou usando o termo técnico, neste último caso, esperamos que o modelo encolher Suas estimativas para a média geral significam mais visivelmente do que no primeiro.
Esse tipo de compartilhamento de informações já é muito útil, mas fica melhor. O modelo girino é um interceptações variadas modelo, como McElreath o chama (ou interceptações aleatóriascomo às vezes é – confuso – chamado) – intercepta Referindo -se à maneira como fazemos previsões para entidades (aqui: tanques), sem variáveis preditivas presentes. Portanto, se pudermos reunir informações sobre interceptações, por que não agrupar informações sobre encostas também? Isso nos permitirá, além disso, fazer uso de relacionamentos Entre variáveis aprendidas em diferentes entidades no conjunto de treinamento.
Então, como você deve ter adivinhado agora, inclinações variadas (ou inclinações aleatóriasse você quiser) é o tópico da postagem de hoje. Novamente, levamos um exemplo do livro de McElreath e mostramos como realizar a mesma coisa com tfprobability.
Café, por favor
Ao contrário do gabinete girino, desta vez trabalhamos com dados simulados. Este é os dados que McElreath usa para introduzir o inclinações variadas técnica de modelagem; Ele então continua e aplica -o a um dos conjuntos de dados mais destacados do livro, o pró-social (ou indiferente, antes!) Chimpanzés. Por hoje, permanecemos com os dados simulados por dois motivos: primeiro, o assunto em si é não trivial o suficiente; e segundo, queremos acompanhar cuidadosamente o que nosso modelo faz e se sua saída está suficientemente próxima dos resultados que McElreath obtido de Stan .
Então, o cenário é esse. Os cafés variam de quão populares são. Em um café well-liked, quando você pede café, é provável que espere. Em um café menos well-liked, você provavelmente será servido muito mais rápido. Isso é uma coisa. Segundo, todos os cafés tendem a estar mais lotados de manhã do que à tarde. Assim, de manhã, você esperará mais do que à tarde – isso vale para os cafés populares e menos populares.
Em termos de interceptações e encostas, podemos imaginar as esperanças da manhã como interceptações, e a tarde resultante aguarda como surgindo devido às encostas das linhas que se juntam todas as manhãs e da tarde, respectivamente.
Então, quando parcialmente pool interceptatemos um “intercepto anterior” (ele mesmo restringido por um anterior, é claro) e um conjunto de interceptações específicas do café que variam em torno dele. Quando parcialmente pool encostastemos uma “inclinação anterior”, refletindo o relacionamento geral entre as esperas da manhã e da tarde e um conjunto de encostas específicas do café, refletindo os relacionamentos individuais. Cognitivamente, isso significa que se você nunca esteve no Café Gerbeaud Em Budapeste, mas já esteve em cafés antes, você pode ter uma idéia menos do que não informada sobre quanto tempo vai esperar; Isso também significa que, se você normalmente pega seu café em seu café favorito de esquina de manhã e agora passa por lá à tarde, tem uma idéia aproximada quanto tempo vai demorar (a saber, menos minutos do que de manhã).
Então isso é tudo? Na verdade, não. Em nosso cenário, interceptações e inclinações estão relacionadas. Se, em um café menos well-liked, eu sempre tomo meu café antes de dois minutos se passaram, há pouco espaço para melhorias. Em um café altamente well-liked, se pudesse demorar facilmente dez minutos de manhã, há um bom potencial para diminuir o tempo de espera da tarde. Portanto, na minha previsão para o tempo de espera desta tarde, devo levar em consideração esse efeito de interação.
Então, agora que temos uma idéia do que se trata, vamos ver como podemos modelar esses efeitos com tfprobability. Mas primeiro, na verdade, temos que gerar os dados.
Simular os dados
Seguimos diretamente McElreath na maneira como os dados são gerados.
##### Inputs wanted to generate the covariance matrix between intercepts and slopes #####
# common morning wait time
a <- 3.5
# common distinction afternoon wait time
# we wait much less within the afternoons
b <- -1
# normal deviation within the (café-specific) intercepts
sigma_a <- 1
# normal deviation within the (café-specific) slopes
sigma_b <- 0.5
# correlation between intercepts and slopes
# the upper the intercept, the extra the wait goes down
rho <- -0.7
##### Generate the covariance matrix #####
# technique of intercepts and slopes
mu <- c(a, b)
# normal deviations of means and slopes
sigmas <- c(sigma_a, sigma_b)
# correlation matrix
# a correlation matrix has ones on the diagonal and the correlation within the off-diagonals
rho <- matrix(c(1, rho, rho, 1), nrow = 2)
# now matrix multiply to get covariance matrix
cov_matrix <- diag(sigmas) %*% rho %*% diag(sigmas)
##### Generate the café-specific intercepts and slopes #####
# 20 cafés total
n_cafes <- 20
library(MASS)
set.seed(5) # used to copy instance
# multivariate distribution of intercepts and slopes
vary_effects <- mvrnorm(n_cafes , mu ,cov_matrix)
# intercepts are within the first column
a_cafe <- vary_effects( ,1)
# slopes are within the second
b_cafe <- vary_effects( ,2)
##### Generate the precise wait occasions #####
set.seed(22)
# 10 visits per café
n_visits <- 10
# alternate values for mornings and afternoons within the knowledge body
afternoon <- rep(0:1, n_visits * n_cafes/2)
# knowledge for every café are consecutive rows within the knowledge body
cafe_id <- rep(1:n_cafes, every = n_visits)
# the regression equation for the imply ready time
mu <- a_cafe(cafe_id) + b_cafe(cafe_id) * afternoon
# normal deviation of ready time inside cafés
sigma <- 0.5 # std dev inside cafes
# generate cases of ready occasions
wait <- rnorm(n_visits * n_cafes, mu, sigma)
d <- knowledge.body(cafe = cafe_id, afternoon = afternoon, wait = wait)Dê um vislumbre dos dados:
Observations: 200
Variables: 3
$ cafe 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,...
$ afternoon 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,...
$ wait 3.9678929, 3.8571978, 4.7278755, 2.7610133, 4.1194827, 3.54365,... Para construir o modelo.
O modelo
Como no Postagem anterior sobre modelagem multinívelnós usamos TFD_JONT_DISTRIBUTION_SEMENTE Para definir o modelo e Hamiltoniano Monte Carlo para amostragem. Considere dar uma olhada na primeira seção dessa postagem para obter um rápido lembrete do procedimento geral.
Antes de codificarmos o modelo, vamos rapidamente tirar a biblioteca do caminho. Importante, novamente como no put up anterior, precisamos instalar um grasp Construição da probabilidade do tensorflow, pois estamos usando recursos muito novos ainda não disponíveis na versão de lançamento atual. O mesmo vale para os pacotes R tensorflow e tfprobability: Por favor, instale as respectivas versões de desenvolvimento do Github.
Agora aqui está a definição do modelo. Vamos passar por isso passo a passo em um instante.
mannequin <- perform(cafe_id) {
tfd_joint_distribution_sequential(
listing(
# rho, the prior for the correlation matrix between intercepts and slopes
tfd_cholesky_lkj(2, 2),
# sigma, prior variance for the ready time
tfd_sample_distribution(tfd_exponential(fee = 1), sample_shape = 1),
# sigma_cafe, prior of variances for intercepts and slopes (vector of two)
tfd_sample_distribution(tfd_exponential(fee = 1), sample_shape = 2),
# b, the prior imply for the slopes
tfd_sample_distribution(tfd_normal(loc = -1, scale = 0.5), sample_shape = 1),
# a, the prior imply for the intercepts
tfd_sample_distribution(tfd_normal(loc = 5, scale = 2), sample_shape = 1),
# mvn, multivariate distribution of intercepts and slopes
# form: batch measurement, 20, 2
perform(a,b,sigma_cafe,sigma,chol_rho)
tfd_sample_distribution(
tfd_multivariate_normal_tri_l(
loc = tf$concat(listing(a,b), axis = -1L),
scale_tril = tf$linalg$LinearOperatorDiag(sigma_cafe)$matmul(chol_rho)),
sample_shape = n_cafes),
# ready time
# form needs to be batch measurement, 200
perform(mvn, a, b, sigma_cafe, sigma)
tfd_independent(
# want to drag out the proper cafe_id within the center column
tfd_normal(
loc = (tf$collect(mvn( , , 1), cafe_id, axis = -1L) +
tf$collect(mvn( , , 2), cafe_id, axis = -1L) * afternoon),
scale=sigma), # Form (batch, 1)
reinterpreted_batch_ndims=1
)
)
)
}As cinco primeiras distribuições são anteriores. Primeiro, temos o anterior para a matriz de correlação. Basicamente, isso seria um Distribuição LKJ de forma 2x2 e com concentração Parâmetro igual a 2.
Por razões de desempenho, trabalhamos com uma versão que insere e produz fatores de Cholesky:
# rho, the prior correlation matrix between intercepts and slopes
tfd_cholesky_lkj(2, 2)Que tipo de antes é esse? Enquanto McElreath continua nos lembrando, nada é mais instrutivo do que amostragem do anterior. Para vermos o que está acontecendo, usamos a distribuição Base LKJ, não a Cholesky:
corr_prior <- tfd_lkj(2, 2)
correlation <- (corr_prior %>% tfd_sample(100))( , 1, 2) %>% as.numeric()
library(ggplot2)
knowledge.body(correlation) %>% ggplot(aes(x = correlation)) + geom_density()
Portanto, esse anterior é moderadamente cético em relação às fortes correlações, mas bastante aberto a aprender com os dados.
A próxima distribuição na fila
# sigma, prior variance for the ready time
tfd_sample_distribution(tfd_exponential(fee = 1), sample_shape = 1)é o anterior para a variação do tempo de espera, a última distribuição na lista.
Em seguida, é a distribuição anterior de variações para as interceptações e inclinações. Este anterior é o mesmo para os dois casos, mas especificamos um sample_shape de 2 para obter duas amostras individuais.
# sigma_cafe, prior of variances for intercepts and slopes (vector of two)
tfd_sample_distribution(tfd_exponential(fee = 1), sample_shape = 2)Agora que temos as respectivas variações anteriores, passamos para os meios anteriores. Ambos são distribuições normais.
# b, the prior imply for the slopes
tfd_sample_distribution(tfd_normal(loc = -1, scale = 0.5), sample_shape = 1)# a, the prior imply for the intercepts
tfd_sample_distribution(tfd_normal(loc = 5, scale = 2), sample_shape = 1)No coração do modelo, onde acontece o pool parcial. Nós vamos construir interceptações e inclinações parcialmente com alças para todos os cafés. Como dissemos acima, interceptações e inclinações não são independentes; Eles interagem. Assim, precisamos usar uma distribuição regular multivariada. Os meios são dados pelos meios anteriores definidos brand acima, enquanto a matriz de covariância é construída a partir das variações anteriores acima e da matriz de correlação anterior. A forma de saída aqui é determinada pelo número de cafés: queremos uma interceptação e uma inclinação para cada café.
# mvn, multivariate distribution of intercepts and slopes
# form: batch measurement, 20, 2
perform(a,b,sigma_cafe,sigma,chol_rho)
tfd_sample_distribution(
tfd_multivariate_normal_tri_l(
loc = tf$concat(listing(a,b), axis = -1L),
scale_tril = tf$linalg$LinearOperatorDiag(sigma_cafe)$matmul(chol_rho)),
sample_shape = n_cafes)Finalmente, provemos os tempos de espera reais. Esse código retira as interceptações e inclinações corretas do regular multivariado e produz o tempo médio de espera, dependendo do café em que estamos e se é de manhã ou da tarde.
# ready time
# form: batch measurement, 200
perform(mvn, a, b, sigma_cafe, sigma)
tfd_independent(
# want to drag out the proper cafe_id within the center column
tfd_normal(
loc = (tf$collect(mvn( , , 1), cafe_id, axis = -1L) +
tf$collect(mvn( , , 2), cafe_id, axis = -1L) * afternoon),
scale=sigma),
reinterpreted_batch_ndims=1
)Antes de executar a amostragem, é sempre uma boa ideia fazer uma verificação rápida no modelo.
n_cafes <- 20
cafe_id <- tf$forged((d$cafe - 1) %% 20, tf$int64)
afternoon <- d$afternoon
wait <- d$waitAmostramos do modelo e, em seguida, verificamos a probabilidade de log.
m <- mannequin(cafe_id)
s <- m %>% tfd_sample(3)
m %>% tfd_log_prob(s)Queremos uma probabilidade de log escalar por membro no lote, que é o que obtemos.
tf.Tensor((-466.1392 -149.92587 -196.51688), form=(3,), dtype=float32)Correndo as correntes
A amostragem actual de Monte Carlo funciona como no put up anterior, com uma exceção. A amostragem acontece no espaço de parâmetros irrestritos, mas no last precisamos obter parâmetros de matriz de correlação válidos rho e variações válidas sigma e sigma_cafe. A conversão entre os espaços é feita por meio de bijectos de TFP. Felizmente, isso não é algo que precisamos fazer como usuários; Tudo o que precisamos especificar são bijetores apropriados. Para as distribuições normais no modelo, não há nada a fazer.
constraining_bijectors <- listing(
# ensure the rho(1:4) parameters are legitimate for a Cholesky issue
tfb_correlation_cholesky(),
# ensure variance is optimistic
tfb_exp(),
# ensure variance is optimistic
tfb_exp(),
tfb_identity(),
tfb_identity(),
tfb_identity()
)Agora podemos montar o amostrador Hamiltonian Monte Carlo.
n_steps <- 500
n_burnin <- 500
n_chains <- 4
# arrange the optimization goal
logprob <- perform(rho, sigma, sigma_cafe, b, a, mvn)
m %>% tfd_log_prob(listing(rho, sigma, sigma_cafe, b, a, mvn, wait))
# preliminary states for the sampling process
c(initial_rho, initial_sigma, initial_sigma_cafe, initial_b, initial_a, initial_mvn, .) %<-%
(m %>% tfd_sample(n_chains))
# HMC sampler, with the above bijectors and step measurement adaptation
hmc <- mcmc_hamiltonian_monte_carlo(
target_log_prob_fn = logprob,
num_leapfrog_steps = 3,
step_size = listing(0.1, 0.1, 0.1, 0.1, 0.1, 0.1)
) %>%
mcmc_transformed_transition_kernel(bijector = constraining_bijectors) %>%
mcmc_simple_step_size_adaptation(target_accept_prob = 0.8,
num_adaptation_steps = n_burnin)Novamente, podemos obter diagnósticos adicionais (aqui: tamanhos de etapas e taxas de aceitação) registrando uma função de rastreamento:
trace_fn <- perform(state, pkr) {
listing(pkr$inner_results$inner_results$is_accepted,
pkr$inner_results$inner_results$accepted_results$step_size)
}Aqui, então, está a função de amostragem. Observe como usamos tf_function para colocá -lo no gráfico. Pelo menos a partir de hoje, isso faz uma enorme diferença no desempenho da amostragem ao usar a execução ansiosa.
run_mcmc <- perform(kernel) {
kernel %>% mcmc_sample_chain(
num_results = n_steps,
num_burnin_steps = n_burnin,
current_state = listing(initial_rho,
tf$ones_like(initial_sigma),
tf$ones_like(initial_sigma_cafe),
initial_b,
initial_a,
initial_mvn),
trace_fn = trace_fn
)
}
run_mcmc <- tf_function(run_mcmc)
res <- hmc %>% run_mcmc()
mcmc_trace <- res$all_statesEntão, como são nossas amostras e o que obtemos em termos de posteriors? Vamos ver.
Resultados
Neste momento, mcmc_trace é uma lista de tensores de formas diferentes, dependentes de como definimos os parâmetros. Precisamos fazer um pouco de pós-processamento para poder resumir e exibir os resultados.
# the precise mcmc samples
# for the hint plots, we need to have them in form (500, 4, 49)
# that's: (variety of steps, variety of chains, variety of parameters)
samples <- abind(
# rho 1:4
as.array(mcmc_trace((1)) %>% tf$reshape(listing(tf$forged(n_steps, tf$int32), tf$forged(n_chains, tf$int32), 4L))),
# sigma
as.array(mcmc_trace((2))),
# sigma_cafe 1:2
as.array(mcmc_trace((3))( , , 1)),
as.array(mcmc_trace((3))( , , 2)),
# b
as.array(mcmc_trace((4))),
# a
as.array(mcmc_trace((5))),
# mvn 10:49
as.array( mcmc_trace((6)) %>% tf$reshape(listing(tf$forged(n_steps, tf$int32), tf$forged(n_chains, tf$int32), 40L))),
alongside = 3)
# the efficient pattern sizes
# we wish them in form (4, 49), which is (variety of chains * variety of parameters)
ess <- mcmc_effective_sample_size(mcmc_trace)
ess <- cbind(
# rho 1:4
as.matrix(ess((1)) %>% tf$reshape(listing(tf$forged(n_chains, tf$int32), 4L))),
# sigma
as.matrix(ess((2))),
# sigma_cafe 1:2
as.matrix(ess((3))( , 1, drop = FALSE)),
as.matrix(ess((3))( , 2, drop = FALSE)),
# b
as.matrix(ess((4))),
# a
as.matrix(ess((5))),
# mvn 10:49
as.matrix(ess((6)) %>% tf$reshape(listing(tf$forged(n_chains, tf$int32), 40L)))
)
# the rhat values
# we wish them in form (49), which is (variety of parameters)
rhat <- mcmc_potential_scale_reduction(mcmc_trace)
rhat <- c(
# rho 1:4
as.double(rhat((1)) %>% tf$reshape(listing(4L))),
# sigma
as.double(rhat((2))),
# sigma_cafe 1:2
as.double(rhat((3))(1)),
as.double(rhat((3))(2)),
# b
as.double(rhat((4))),
# a
as.double(rhat((5))),
# mvn 10:49
as.double(rhat((6)) %>% tf$reshape(listing(40L)))
) Lotes de rastreamento
Quão bem as correntes se misturam?
prep_tibble <- perform(samples) {
as_tibble(samples, .name_repair = ~ c("chain_1", "chain_2", "chain_3", "chain_4")) %>%
add_column(pattern = 1:n_steps) %>%
collect(key = "chain", worth = "worth", -pattern)
}
plot_trace <- perform(samples) {
prep_tibble(samples) %>%
ggplot(aes(x = pattern, y = worth, shade = chain)) +
geom_line() +
theme_light() +
theme(legend.place = "none",
axis.title = element_blank(),
axis.textual content = element_blank(),
axis.ticks = element_blank())
}
plot_traces <- perform(sample_array, num_params) {
plots <- purrr::map(1:num_params, ~ plot_trace(sample_array( , , .x)))
do.name(grid.prepare, plots)
}
plot_traces(samples, 49)
Incrível! (Os dois primeiros parâmetros de rhoo fator Cholesky da matriz de correlação, precisa permanecer fixo em 1 e 0, respectivamente.)
Agora, para algumas estatísticas resumidas sobre os posteriors dos parâmetros.
Parâmetros
Como na última vez, exibimos meios posteriores e desvios padrão, bem como o maior intervalo de densidade posterior (HPDI). Adicionamos tamanhos de amostra eficazes e rhat valores.
column_names <- c(
paste0("rho_", 1:4),
"sigma",
paste0("sigma_cafe_", 1:2),
"b",
"a",
c(rbind(paste0("a_cafe_", 1:20), paste0("b_cafe_", 1:20)))
)
all_samples <- matrix(samples, nrow = n_steps * n_chains, ncol = 49)
all_samples <- all_samples %>%
as_tibble(.name_repair = ~ column_names)
all_samples %>% glimpse()
means <- all_samples %>%
summarise_all(listing (imply)) %>%
collect(key = "key", worth = "imply")
sds <- all_samples %>%
summarise_all(listing (sd)) %>%
collect(key = "key", worth = "sd")
hpdis <-
all_samples %>%
summarise_all(listing(~ listing(hdi(.) %>% t() %>% as_tibble()))) %>%
unnest()
hpdis_lower <- hpdis %>% choose(-accommodates("higher")) %>%
rename(lower0 = decrease) %>%
collect(key = "key", worth = "decrease") %>%
prepare(as.integer(str_sub(key, 6))) %>%
mutate(key = column_names)
hpdis_upper <- hpdis %>% choose(-accommodates("decrease")) %>%
rename(upper0 = higher) %>%
collect(key = "key", worth = "higher") %>%
prepare(as.integer(str_sub(key, 6))) %>%
mutate(key = column_names)
abstract <- means %>%
inner_join(sds, by = "key") %>%
inner_join(hpdis_lower, by = "key") %>%
inner_join(hpdis_upper, by = "key")
ess <- apply(ess, 2, imply)
summary_with_diag <- abstract %>% add_column(ess = ess, rhat = rhat)
print(summary_with_diag, n = 49)# A tibble: 49 x 7
key imply sd decrease higher ess rhat
1 rho_1 1 0 1 1 NaN NaN
2 rho_2 0 0 0 0 NaN NaN
3 rho_3 -0.517 0.176 -0.831 -0.195 42.4 1.01
4 rho_4 0.832 0.103 0.644 1.000 46.5 1.02
5 sigma 0.473 0.0264 0.420 0.523 424. 1.00
6 sigma_cafe_1 0.967 0.163 0.694 1.29 97.9 1.00
7 sigma_cafe_2 0.607 0.129 0.386 0.861 42.3 1.03
8 b -1.14 0.141 -1.43 -0.864 95.1 1.00
9 a 3.66 0.218 3.22 4.07 75.3 1.01
10 a_cafe_1 4.20 0.192 3.83 4.57 83.9 1.01
11 b_cafe_1 -1.13 0.251 -1.63 -0.664 63.6 1.02
12 a_cafe_2 2.17 0.195 1.79 2.54 59.3 1.01
13 b_cafe_2 -0.923 0.260 -1.42 -0.388 46.0 1.01
14 a_cafe_3 4.40 0.195 4.02 4.79 56.7 1.01
15 b_cafe_3 -1.97 0.258 -2.52 -1.51 43.9 1.01
16 a_cafe_4 3.22 0.199 2.80 3.57 58.7 1.02
17 b_cafe_4 -1.20 0.254 -1.70 -0.713 36.3 1.01
18 a_cafe_5 1.86 0.197 1.45 2.20 52.8 1.03
19 b_cafe_5 -0.113 0.263 -0.615 0.390 34.6 1.04
20 a_cafe_6 4.26 0.210 3.87 4.67 43.4 1.02
21 b_cafe_6 -1.30 0.277 -1.80 -0.713 41.4 1.05
22 a_cafe_7 3.61 0.198 3.23 3.98 44.9 1.01
23 b_cafe_7 -1.02 0.263 -1.51 -0.489 37.7 1.03
24 a_cafe_8 3.95 0.189 3.59 4.31 73.1 1.01
25 b_cafe_8 -1.64 0.248 -2.10 -1.13 60.7 1.02
26 a_cafe_9 3.98 0.212 3.57 4.37 76.3 1.03
27 b_cafe_9 -1.29 0.273 -1.83 -0.776 57.8 1.05
28 a_cafe_10 3.60 0.187 3.24 3.96 104. 1.01
29 b_cafe_10 -1.00 0.245 -1.47 -0.512 70.4 1.00
30 a_cafe_11 1.95 0.200 1.56 2.35 55.9 1.03
31 b_cafe_11 -0.449 0.266 -1.00 0.0619 42.5 1.04
32 a_cafe_12 3.84 0.195 3.46 4.22 76.0 1.02
33 b_cafe_12 -1.17 0.259 -1.65 -0.670 62.5 1.03
34 a_cafe_13 3.88 0.201 3.50 4.29 62.2 1.02
35 b_cafe_13 -1.81 0.270 -2.30 -1.29 48.3 1.03
36 a_cafe_14 3.19 0.212 2.82 3.61 65.9 1.07
37 b_cafe_14 -0.961 0.278 -1.49 -0.401 49.9 1.06
38 a_cafe_15 4.46 0.212 4.08 4.91 62.0 1.09
39 b_cafe_15 -2.20 0.290 -2.72 -1.59 47.8 1.11
40 a_cafe_16 3.41 0.193 3.02 3.78 62.7 1.02
41 b_cafe_16 -1.07 0.253 -1.54 -0.567 48.5 1.05
42 a_cafe_17 4.22 0.201 3.82 4.60 58.7 1.01
43 b_cafe_17 -1.24 0.273 -1.74 -0.703 43.8 1.01
44 a_cafe_18 5.77 0.210 5.34 6.18 66.0 1.02
45 b_cafe_18 -1.05 0.284 -1.61 -0.511 49.8 1.02
46 a_cafe_19 3.23 0.203 2.88 3.65 52.7 1.02
47 b_cafe_19 -0.232 0.276 -0.808 0.243 45.2 1.01
48 a_cafe_20 3.74 0.212 3.35 4.21 48.2 1.04
49 b_cafe_20 -1.09 0.281 -1.58 -0.506 36.5 1.05 Então, o que temos? Se você executar este “ao vivo”, para as linhas a_cafe_n resp. b_cafe_nvocê vê uma boa alternância de coloração branca e vermelha: para todos os cafés, as encostas inferidas são negativas.
A inclinação inferida antes (b) está em torno de -1,14, que não está muito longe do valor que usamos para amostragem: 1.
O rho As estimativas posteriores, reconhecidamente, são menos úteis, a menos que você esteja acostumado a compor fatores de Cholesky em sua cabeça. Nós calculamos as correlações posteriores resultantes e sua média:
-0.5166775O valor que usamos para amostragem foi -0,7, por isso vemos o efeito de regularização. Caso você esteja se perguntando, pelos mesmos dados Stan produz uma estimativa de -0,5.
Finalmente, vamos exibir equivalentes às figuras de McElreath que ilustram o encolhimento no parâmetro (interceptações e inclinações específicas do café), bem como o resultado (tempos de espera da tarde da tarde).
Encolhimento
Como esperado, vemos que as interceptações e inclinações individuais são puxadas para a média – quanto mais, mais longe são do centro.
# identical to McElreath, compute unpooled estimates straight from knowledge
a_empirical <- d %>%
filter(afternoon == 0) %>%
group_by(cafe) %>%
summarise(a = imply(wait)) %>%
choose(a)
b_empirical <- d %>%
filter(afternoon == 1) %>%
group_by(cafe) %>%
summarise(b = imply(wait)) %>%
choose(b) -
a_empirical
empirical_estimates <- bind_cols(
a_empirical,
b_empirical,
kind = rep("knowledge", 20))
posterior_estimates <- tibble(
a = means %>% filter(
str_detect(key, "^a_cafe")) %>% choose(imply) %>% pull(),
b = means %>% filter(
str_detect(key, "^b_cafe")) %>% choose(imply) %>% pull(),
kind = rep("posterior", 20))
all_estimates <- bind_rows(empirical_estimates, posterior_estimates)
# compute posterior imply bivariate Gaussian
# once more following McElreath
mu_est <- c(means(means$key == "a", 2), means(means$key == "b", 2)) %>% unlist()
rho_est <- mean_rho
sa_est <- means(means$key == "sigma_cafe_1", 2) %>% unlist()
sb_est <- means(means$key == "sigma_cafe_2", 2) %>% unlist()
cov_ab <- sa_est * sb_est * rho_est
sigma_est <- matrix(c(sa_est^2, cov_ab, cov_ab, sb_est^2), ncol=2)
alpha_levels <- c(0.1, 0.3, 0.5, 0.8, 0.99)
names(alpha_levels) <- alpha_levels
contour_data <- plyr::ldply(
alpha_levels,
ellipse,
x = sigma_est,
scale = c(1, 1),
centre = mu_est
)
ggplot() +
geom_point(knowledge = all_estimates, mapping = aes(x = a, y = b, shade = kind)) +
geom_path(knowledge = contour_data, mapping = aes(x = x, y = y, group = .id))
O mesmo comportamento é visível na escala de resultados.
wait_times <- all_estimates %>%
mutate(morning = a, afternoon = a + b)
# simulate from posterior means
v <- MASS::mvrnorm(1e4 , mu_est , sigma_est)
v( ,2) <- v( ,1) + v( ,2) # calculate afternoon wait
# assemble empirical covariance matrix
sigma_est2 <- cov(v)
mu_est2 <- mu_est
mu_est2(2) <- mu_est(1) + mu_est(2)
contour_data <- plyr::ldply(
alpha_levels,
ellipse,
x = sigma_est2 %>% unname(),
scale = c(1, 1),
centre = mu_est2
)
ggplot() +
geom_point(knowledge = wait_times, mapping = aes(x = morning, y = afternoon, shade = kind)) +
geom_path(knowledge = contour_data, mapping = aes(x = x, y = y, group = .id))
Embrulhando
Até agora, esperamos que tenhamos convencido você do poder inerente à modelagem bayesiana, além de transmitir algumas idéias sobre como isso é possível com a probabilidade de tensorflow. Como em cada DSL, porém, leva tempo para proceder da compreensão de exemplos trabalhados para projetar seus próprios modelos. E não apenas tempo – ajuda a ter visto muitos modelos diferentes, concentrando -se em diferentes tarefas e aplicativos. Neste weblog, planejamos acompanhar vagamente a modelagem bayesiana com o TFP, pegando algumas das tarefas e desafios elaborados nos capítulos posteriores do livro de McElreath. Obrigado pela leitura!