Estudo Simulado do modelo de Regressão Dirichlet

Abordagem Bayesiana - Parte 5

Pedro Frazão Dutra

20/07/2026

O presente documento apresenta os resultados da aplicação do modelo de Regressão Dirichlet a dados simulados. O cenário aqui considerado foi uma amostra com tamanho \(n = 750\) de vetores composicionais de dimensão 7. As covariáveis associadas possuem, respectivamente, distribuição uniforme no conjunto \([0,10]\) e distribuição Bernoulli com parâmetro \(p = 1/2\), ou seja, geramos:

\[Y_i = (Y_{i1}, Y_{i2}, Y_{i3}, Y_{i4}, Y_{i5}, Y_{i6}, Y_{i7}) \sim \text{Dirichlet}(\alpha_{i1}, \alpha_{i2}, \alpha_{i3}, \alpha_{i4}, \alpha_{i5}, \alpha_{i6}, \alpha_{i7})\]

\[\alpha_{ij} = \exp(\beta_{0j} + \beta_{1j} X_{1i} + \beta_{2j} X_{2i}), \quad \text{para } j = 1, \dots, 7\]

\[\beta \sim N(0, 100 I_d)\]

Para a simulação, os valores verdadeiros dos parâmetros foram organizados na matriz de coeficientes \(\boldsymbol{\beta}\), onde cada linha representa o efeito de uma covariável e cada coluna mapeia uma das componentes da resposta composicional:

\[\boldsymbol{\beta} = \begin{bmatrix} \beta_{01} & \beta_{02} & \beta_{03} & \beta_{04} & \beta_{05} & \beta_{06} & \beta_{07} \\ \beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} & \beta_{15} & \beta_{16} & \beta_{17} \\ \beta_{21} & \beta_{22} & \beta_{23} & \beta_{24} & \beta_{25} & \beta_{26} & \beta_{27} \end{bmatrix} = \begin{bmatrix} \phantom{-}0.5 & \phantom{-}0.0 & -0.5 & \phantom{-}0.2 & -0.4 & \phantom{-}0.1 & -0.2 \\ -0.8 & \phantom{-}0.6 & \phantom{-}0.3 & -0.1 & \phantom{-}0.5 & -0.4 & \phantom{-}0.2 \\ \phantom{-}1.2 & -0.9 & \phantom{-}0.4 & -0.7 & \phantom{-}0.2 & \phantom{-}0.8 & -0.3 \end{bmatrix}\]

Onde: * Linha 1 (\(\beta_{0j}\)): Interceptos para as componentes \(Y_1, Y_2, Y_3, Y_4, Y_5, Y_6\) e \(Y_7\). * Linha 2 (\(\beta_{1j}\)): Efeitos da covariável \(X_1\) sobre as componentes. * Linha 3 (\(\beta_{2j}\)): Efeitos da covariável \(X_2\) sobre as componentes.

1 Pacotes utilizados

# Análise composicional
library(DirichletReg) 
library(compositions) 

# Produção de gráficos e tabelas
library(tidyverse)
library(ggtern)
library(colorspace)
library(purrr)
library(patchwork)
library(DT)          
library(knitr)       
library(kableExtra)  
library(htmltools)   
library(scales)

# Inferência Bayesiana
library(posterior)
library(coda)

2 Geração dos dados

set.seed(2005)

# Tamanho da amostra
n <- 750

# Simulando e centralizando a covariável X1 (Uniforme)
x1_bruto <- runif(n, 0, 5)
x1 <- x1_bruto - mean(x1_bruto)

# Simulando a covariável X2 (Bernoulli com p = 0.5)
x2 <- rbinom(n, 1, 0.5)

# Matriz de design X (Intercepto, x1 e x2)
X_mat <- cbind(1, x1, x2)

# Matriz de Betas 
beta_real <- matrix(c(
   0.5,  0.0, -0.5,  0.2, -0.4,  0.1, -0.2,  # Interceptos
  -0.8,  0.6,  0.3, -0.1,  0.5, -0.4,  0.2,  # Efeitos de X1 
   1.2, -0.9,  0.4, -0.7,  0.2,  0.8, -0.3   # Efeitos de X2 
), nrow = 3, byrow = TRUE)

# Função de ligação
log_alpha <- X_mat %*% beta_real

# Cálculo dos alphas
alpha <- exp(log_alpha)

# Gerando as composições 
y <- matrix(0, nrow = n, ncol = 7)

for(i in 1:n) {
  z <- rgamma(7, shape = alpha[i, ], rate = 1)
  
  # Fecha no simplex (soma 1)
  y[i, ] <- z / sum(z)
}

# Expandindo os nomes dos parâmetros para as 7 dimensões 
nomes_param <- c("y1: Intercepto", "y1: Slope v1", "y1: Slope v2",
                 "y2: Intercepto", "y2: Slope v1", "y2: Slope v2",
                 "y3: Intercepto", "y3: Slope v1", "y3: Slope v2",
                 "y4: Intercepto", "y4: Slope v1", "y4: Slope v2",
                 "y5: Intercepto", "y5: Slope v1", "y5: Slope v2",
                 "y6: Intercepto", "y6: Slope v1", "y6: Slope v2",
                 "y7: Intercepto", "y7: Slope v1", "y7: Slope v2")

3 Análise descritiva

4 Ajuste clássico do modelo do Maier

## Call:
## DirichReg(formula = AL ~ x1 + x2, data = dados_modelo)
## 
## Standardized Residuals:
##         Min       1Q   Median      3Q      Max
## v1  -3.7902  -0.6920  -0.0673  0.6947   4.8707
## v2  -1.8957  -0.6359  -0.3969  0.2460  11.2788
## v3  -1.4190  -0.7009  -0.3606  0.3347   5.7084
## v4  -1.3720  -0.7414  -0.3107  0.4664   4.0343
## v5  -1.6185  -0.6764  -0.3436  0.4178   5.7022
## v6  -1.9590  -0.7791  -0.2396  0.6212   4.7815
## v7  -1.1837  -0.7156  -0.3760  0.4252   5.4053
## 
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 1: v1
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.58432    0.04256   13.73   <2e-16 ***
## x1          -0.82610    0.01985  -41.61   <2e-16 ***
## x2           1.21104    0.05547   21.83   <2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 2: v2
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.09084    0.04502   2.018   0.0436 *  
## x1           0.49454    0.02332  21.208   <2e-16 ***
## x2          -0.72734    0.06400 -11.364   <2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 3: v3
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.43769    0.04853  -9.020  < 2e-16 ***
## x1           0.24182    0.02249  10.753  < 2e-16 ***
## x2           0.41839    0.06525   6.413 1.43e-10 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 4: v4
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.32873    0.04456   7.377 1.61e-13 ***
## x1          -0.13734    0.02231  -6.155 7.50e-10 ***
## x2          -0.74122    0.06388 -11.604  < 2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 5: v5
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.26232    0.04716  -5.562 2.67e-08 ***
## x1           0.40142    0.02267  17.710  < 2e-16 ***
## x2           0.21569    0.06386   3.378 0.000731 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 6: v6
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.19534    0.04483   4.357 1.32e-05 ***
## x1          -0.42290    0.02052 -20.613  < 2e-16 ***
## x2           0.79014    0.05872  13.455  < 2e-16 ***
## ------------------------------------------------------------------
## Beta-Coefficients for variable no. 7: v7
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.12569    0.04708  -2.670 0.007586 ** 
## x1           0.15309    0.02286   6.696 2.14e-11 ***
## x2          -0.24819    0.06556  -3.786 0.000153 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 8210 on 21 df (138 BFGS + 2 NR Iterations)
## AIC: -16378, BIC: -16281
## Number of Observations: 750
## Link: Log
## Parametrization: common

5 Rotina de Metropolis-Hastings

5.1 Funções de Log-Verossimilhança e Log-Posteriori

log_verossimilhanca <- function(beta_vec, Y, X) {
  N <- nrow(Y)
  C <- ncol(Y)
  K <- ncol(X)
  
  beta_mat <- matrix(beta_vec, nrow = K, ncol = C)
  alpha <- exp(X %*% beta_mat)
  
  termo1 <- sum(lgamma(rowSums(alpha)))
  termo2 <- sum(lgamma(alpha))
  termo3 <- sum((alpha - 1) * log(Y))
  
  ll <- termo1 - termo2 + termo3
  return(ll)
}

log_posteriori <- function(beta_vec, Y, X, sigma_prior = 100) {
  ll <- log_verossimilhanca(beta_vec, Y, X)
  log_prior <- sum(dnorm(beta_vec, mean = 0, sd = sigma_prior, log = TRUE))
  log_post <- ll + log_prior
  return(log_post)
}

5.2 Algoritmo

dados_modelo <- data.frame(
  x1 = x1, 
  x2 = x2, 
  y1 = y[, 1], y2 = y[, 2], y3 = y[, 3], 
  y4 = y[, 4], y5 = y[, 5], y6 = y[, 6], y7 = y[, 7]
)

C <- 7 
K <- 3 
d <- C * K # 21 dimensões

Y_matriz <- as.matrix(dados_modelo[, c("y1","y2","y3","y4","y5","y6","y7")])
X_matriz <- cbind(1, dados_modelo$x1, dados_modelo$x2)

# Configurações para múltiplas cadeias
m <- 4          
nite <- 100000  

# Parâmetros de base do modelo clássico
beta_mle <- as.numeric(unlist(coef(modelo)))
Sigma_mle <- vcov(modelo)

# Fator de tuning recalibrado para a geometria de 21 dimensões
tuning <- 1 
Sigma_prop <- tuning * (2.4^2 / d) * Sigma_mle

# Inicializando lista para armazenar o histórico de cada cadeia
cadeias <- list()
set.seed(1005) 

# Algoritmo Otimizado para Alta Dimensão
for (i in 1:m) {
  
  tempo_inicio <- Sys.time()
  
  # Inicialização controlada: evita o "vazio" de densidade em 21D
  beta_ini <- MASS::mvrnorm(1, mu = beta_mle, Sigma = 5 * Sigma_mle)
  
  beta_cadeia <- matrix(0, nrow = nite, ncol = d)
  beta_cadeia[1, ] <- beta_ini
  
  log_post_atual <- log_posteriori(beta_ini, Y = Y_matriz, X = X_matriz)
  
  # Pré-simulação dos ruídos (vetorização para ganho de performance)
  ruido_prop <- MASS::mvrnorm(nite - 1, mu = rep(0, d), Sigma = Sigma_prop)
  
  # Passo de Metropolis
  aceitos <- 0
  for (t in 1:(nite - 1)) {
    
    beta_prop <- beta_cadeia[t, ] + ruido_prop[t, ]
    
    log_post_prop <- log_posteriori(beta_prop, Y = Y_matriz, X = X_matriz)
    
    # Proteção robusta contra indefinições numéricas num espaço 21-dimensional
    if (is.na(log_post_prop) || is.infinite(log_post_prop)) {
      log_alfa <- -Inf 
    } else {
      log_alfa <- log_post_prop - log_post_atual
    }
    
    if (!is.na(log_alfa) && log(runif(1)) <= log_alfa) {
      beta_cadeia[t+1, ] <- beta_prop
      log_post_atual <- log_post_prop
      aceitos <- aceitos + 1
    } else {
      beta_cadeia[t+1, ] <- beta_cadeia[t, ]
    }
  }
  
  cadeias[[i]] <- beta_cadeia
  
  tempo_fim <- Sys.time()
  tempo_execucao <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)
  
  cat("Cadeia", i, 
      "- Taxa de aceitação:", round((aceitos / (nite - 1)) * 100, 2), "%",
      "| Tempo de execução:", tempo_execucao, "segundos\n")
}
## Cadeia 1 - Taxa de aceitação: 24.87 % | Tempo de execução: 93.3 segundos
## Cadeia 2 - Taxa de aceitação: 24.99 % | Tempo de execução: 90.9 segundos
## Cadeia 3 - Taxa de aceitação: 24.84 % | Tempo de execução: 95.35 segundos
## Cadeia 4 - Taxa de aceitação: 25.06 % | Tempo de execução: 93.11 segundos

5.3 Diagnósticos de convergência

5.3.1 Cadeias piloto

Nesta etapa, simulamos 4 cadeias para cada um dos parâmetos, todas com inicialização superdispersa. Com isso, podemos checar se existe multimodalidade em cada uma das distribuições posteriori. Além disso, para avaliar a convergência de cada grupo de cadeias para a uma mesma distribuição comum, calculamos a estatística potencial de redução de escala, \(\hat{R}\).

5.3.2 Funções de Autocorrelação

Plotamos as funções de Autocorrelação das cadeias para verificar a eficiência com que o algoritmo explora o espaço paramétrico da distribuições posteriori. Além disso, é importante descobrir o menor lag k tal que as ACF’s apresentem todos os seus valores dentro do intervalo de confiança centrado em zero. Com isso, se for preciso, poderemos efetuar um espaçamento de tamanho \(k\) em todas as cadeia, de modo a restarem somente amostras estatísticamente independentes para realizarmos estimações. Como existem muitos parâmetros, plotamos apenas as ACF’s da primeira cadeia. Isto é adequado, uma vez que as quatro cadeias apresentaram um comportamento empíricamente semelhante.

5.3.3 R-hat & ESS

Por fim, calculamos a estatística de redução de escala potencial, \(\hat{R}\) e o tamanho efetivo das amostras para estimação de quantidades próximas do centro de massa da distribuição posteriori e para estimação de quantidades próximas das caudas. Para cada parâmetro, juntamos as cadeias aquecidas e realizamos ambos os cálculos.

Diagnósticos de Convergência
Parâmetro R-Hat ESS Bulk ESS Tail
y1: Intercepto 1.001 3173 6158
y1: Slope v1 1.001 3090 5723
y1: Slope v2 1.001 3147 6645
y2: Intercepto 1.001 2686 5551
y2: Slope v1 1.001 3271 6913
y2: Slope v2 1.001 3064 6080
y3: Intercepto 1.002 3170 5996
y3: Slope v1 1.002 3283 6254
y3: Slope v2 1.002 3090 6159
y4: Intercepto 1.000 2965 6336
y4: Slope v1 1.002 3013 6407
y4: Slope v2 1.001 3044 6517
y5: Intercepto 1.003 3185 6155
y5: Slope v1 1.002 3113 6255
y5: Slope v2 1.001 3421 6675
y6: Intercepto 1.001 2881 5467
y6: Slope v1 1.001 3094 6191
y6: Slope v2 1.002 2758 5240
y7: Intercepto 1.000 3125 6323
y7: Slope v1 1.002 3277 6011
y7: Slope v2 1.000 3057 6100

5.4 Inferência Bayesiana

5.4.1 Distribuições marginais

Assumido com segurança o bom desempenho do algoritmo, plotamos os histogramas das distribuições marginais e calculamos medidas resumos usuais, como média, mediana e intervalos de 95% de credibilidade.

Estimativas Pontuais e Intervalos de 95% de Credibilidade
Parâmetro Valor Real Média Post. Mediana Desvio Padrão Quantil 2.5% Quantil 97.5%
y1: Intercepto 0.5 0.509 0.509 0.043 0.424 0.592
y1: Slope v1 -0.8 -0.804 -0.804 0.020 -0.844 -0.765
y1: Slope v2 1.2 1.198 1.197 0.056 1.090 1.308
y2: Intercepto 0.0 -0.025 -0.025 0.046 -0.116 0.065
y2: Slope v1 0.6 0.616 0.616 0.022 0.572 0.660
y2: Slope v2 -0.9 -0.875 -0.875 0.065 -1.002 -0.749
y3: Intercepto -0.5 -0.508 -0.508 0.048 -0.605 -0.415
y3: Slope v1 0.3 0.292 0.292 0.022 0.248 0.336
y3: Slope v2 0.4 0.398 0.399 0.066 0.269 0.527
y4: Intercepto 0.2 0.289 0.289 0.045 0.201 0.375
y4: Slope v1 -0.1 -0.110 -0.110 0.023 -0.154 -0.067
y4: Slope v2 -0.7 -0.776 -0.776 0.064 -0.900 -0.652
y5: Intercepto -0.4 -0.350 -0.349 0.047 -0.443 -0.259
y5: Slope v1 0.5 0.474 0.474 0.022 0.430 0.517
y5: Slope v2 0.2 0.184 0.184 0.063 0.062 0.308
y6: Intercepto 0.1 0.141 0.142 0.046 0.050 0.229
y6: Slope v1 -0.4 -0.397 -0.396 0.021 -0.437 -0.357
y6: Slope v2 0.8 0.768 0.767 0.060 0.653 0.887
y7: Intercepto -0.2 -0.165 -0.164 0.048 -0.259 -0.071
y7: Slope v1 0.2 0.189 0.189 0.023 0.144 0.235
y7: Slope v2 -0.3 -0.285 -0.285 0.065 -0.416 -0.157

6 Avaliando a qualidade do ajuste do modelo

6.1 Métricas de desempenho

Nesta seção, constam os códigos referente a amostragem da distribuição das 3 métricas de desempenho do modelo utilizadas: Distância de Aitchson, Erro quadrático Médio Padrão e Divergência de Kullback-Leibler.

S <- nrow(cadeia_thin_df)
n <- nrow(y)
K <- 3  # Número de parâmetros
C <- 7  # Número de componentes

# Vetores para guardar as distribuições das métricas
distancias_aitchison <- numeric(S)
rmse_amostras        <- numeric(S)
kl_amostras          <- numeric(S)

# Matrizes para acumular a análise preditiva
Y_hat_acumulado <- matrix(0, nrow = n, ncol = C)
Y_sim_acumulado <- matrix(0, nrow = n, ncol = C)
Y_sim_unica     <- matrix(0, nrow = n, ncol = C) 

matriz_betas_loop <- as.matrix(cadeia_thin_df[, nomes_param])

# ---------------------------------------------------------------------------
# PRÉ-COMPUTAÇÃO (FORA DO LOOP): Transformação CLR dos dados observados Y
# ----------------------------------------------------------------------------
log_y_estavel <- log(y + 1e-10)
clr_y <- log_y_estavel - rowMeans(log_y_estavel)

# CRONÔMETRO: Dispara o relógio antes de entrar no loop
tempo_inicio <- Sys.time()

semente_da_simulacao <- 5555

# --- O SUPER LOOP MESTRE (VERSÃO VETORIZADA - SEM LOOPS INTERNOS)
for (s in 1:S) {
  set.seed(semente_da_simulacao + s)
  
  # PASSO A: Reconstrução dos Parâmetros
  vetor_betas_s <- matriz_betas_loop[s, ]
  beta_s        <- matrix(vetor_betas_s, nrow = K, ncol = C, byrow = FALSE)
  
  # MÉTRICA 1: RMSE (Espaço dos Parâmetros)
  rmse_amostras[s] <- sqrt(mean((beta_real - beta_s)^2))
  
  # PASSO B: Projeção no Simplex
  alpha_s <- exp(X_mat %*% beta_s)
  Y_hat_s <- alpha_s / rowSums(alpha_s)
  
  # Acumula para calcular a tendência média no final
  Y_hat_acumulado <- Y_hat_acumulado + Y_hat_s
  
  # -------------------------------------------------------------------------
  # MÉTRICA 2: Distância de Aitchison Vetorizada 
  # -------------------------------------------------------------------------
  log_Y_hat_estavel <- log(Y_hat_s + 1e-10)
  clr_Y_hat <- log_Y_hat_estavel - rowMeans(log_Y_hat_estavel)
  distancias_aitchison[s] <- mean(sqrt(rowSums((clr_y - clr_Y_hat)^2)))
  
  # -------------------------------------------------------------------------
  # MÉTRICA 3: Divergência Kullback-Leibler
  # -------------------------------------------------------------------------
  kl_individual   <- rowSums(y * log((y + 1e-10) / (Y_hat_s + 1e-10)))
  kl_amostras[s]  <- mean(kl_individual)
  
  # -------------------------------------------------------------------------
  # ANÁLISE PREDITIVA: Amostragem Dirichlet via Gamma Vetorizada 
  # -------------------------------------------------------------------------
  gamas_sim <- matrix(rgamma(n * C, shape = alpha_s, rate = 1), nrow = n, ncol = C)
  Y_sim_s   <- gamas_sim / rowSums(gamas_sim)
  
  # Acumula a simulação e guarda a última para o gráfico de dispersão
  Y_sim_acumulado <- Y_sim_acumulado + Y_sim_s
  if (s == S) {
    Y_sim_unica <- Y_sim_s
  }
}

# CRONÔMETRO: Para o relógio imediatamente após o fim do loop
tempo_fim <- Sys.time()
tempo_loop <- round(difftime(tempo_fim, tempo_inicio, units = "secs"), 2)

# --- EXTRAÇÃO DAS MÉDIAS GLOBAIS ---
distancia_media_aitchison <- mean(distancias_aitchison)
rmse_global               <- mean(rmse_amostras)
kl_media                  <- mean(kl_amostras)

Y_hat_medio <- Y_hat_acumulado / S  
Y_sim_medio <- Y_sim_acumulado / S  

# Exibe o tempo de execução formatado
cat("Tempo de execução do Super Loop:", tempo_loop, "segundos\n")
## Tempo de execução do Super Loop: 2.74 segundos

6.1.1 Distância de Aitchson

6.1.2 Erro quadrático médio padrão

6.2 Divergência de Kullback-Leibler

6.3 Análise preditiva a posteriori

Aqui, simulados dados de acordo com os parâmetros gerados pelo modelo e comparamos com os dados observados

# 1. Garante a nomeação correta para as 7 componentes atuais
colnames(y) <- paste0("Y", 1:7)
colnames(Y_sim_unica) <- paste0("Y", 1:7)

# 2. Estrutura os dados observados e simulados
df_obs <- as.data.frame(y) %>% mutate(Tipo = "Observado")
df_sim <- as.data.frame(Y_sim_unica) %>% mutate(Tipo = "Simulado")

# 3. Unifica e transforma para o formato longo (long format)
df_ppc_long <- rbind(df_obs, df_sim) %>%
  pivot_longer(
    cols = starts_with("Y"),
    names_to = "Componente",
    values_to = "Proporcao"
  )

# 4. Gráfico de Distribuição Facetado (Violin + Boxplot)
grafico_ppc_dist <- ggplot(df_ppc_long, aes(x = Tipo, y = Proporcao, fill = Tipo)) +
  # O violin mostra o formato da densidade dos dados
  geom_violin(alpha = 0.5, color = NA, position = position_dodge(0.8)) +
  # O boxplot traz os quartis e mediana de forma direta
  geom_boxplot(width = 0.22, color = "grey20", outlier.alpha = 0.1, 
               position = position_dodge(0.8), lwd = 0.4) +
  
  # Faceta para cada uma das 7 componentes
  facet_wrap(~ Componente, ncol = 4, scales = "free_y") +
  
  # Cores com alto contraste e fáceis de ler
  scale_fill_manual(values = c("Observado" = "#d95f02", "Simulado" = "#7570b3")) +
  
  theme_minimal(base_size = 11) +
  labs(
    title = "Checagem Preditiva Posterior (PPC)",
    subtitle = "Comparação das distribuições das proporções observadas vs. simuladas",
    x = NULL,
    y = "Proporção (Simplex)",
    fill = "Origem dos Dados"
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
    plot.subtitle = element_text(size = 10, color = "grey40", hjust = 0.5),
    strip.background = element_rect(fill = "gray95", color = NA),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

print(grafico_ppc_dist)