(Acest articol a fost publicat pentru prima dată pe R funcționeazăși a contribuit cu drag la R-Bloggers). (Puteți raporta problema despre conținutul de pe această pagină aici)
Doriți să vă împărtășiți conținutul pe R-Bloggers? Faceți clic aici dacă aveți un blog sau aici dacă nu.
În 1839, matematicianul talentat Peter Gustav Lejeune Dirichlet a fost atașat de Departamentul de Filozofie de la Universitatea din Berlin, care lucra pentru mai puțin decât salariul complet, chiar dacă a devenit membru al Academiei Prusoane de Științe în 1832. Habilitațiischrift Prelegere în latină. Aparent, facilitatea lui Dirichlet cu latină nu a fost în funcție de sarcină, așa că, la fel ca mulți „profesori adjuncți” pricepuți astăzi, Dirichlet a luat un concert lateral pentru a -și susține familia. A predat matematica la o școală militară. Oricum, mă descurc. Aproape de acea vreme Dirichlet a început să lucreze la o problemă în mecanica cerească care a implicat această integrală:
Aici 
care este atras de un punct 
unde 
este forța atracției 
şi 
este norma euclidiană.
După o serie de manipulări din seria supranatural, detaliate de Gupta și Richards, Dirichlet a ajuns la următoarea integrală
pe care o veți recunoaște ca funcție beta, constanta de normalizare pentru distribuția Dirichlet:
cu medie 
și variație 
unde
şi 
Distribuția Dirichlet este o generalizare multivariată a distribuției beta, care este adesea utilizată în statisticile bayesiene ca distribuție prealabilă pentru distribuțiile categorice și multinomiale. Am ilustrat această utilizare a Dirichletului într-un post anterior, în timp ce am construit un model Bayesian pentru un lanț Markov cu trei state. Distribuția Dirichlet este remarcabilă prin faptul că reunește lucrările din secolul al XVIII -lea și al XIX -lea în analiză, așa cum este exemplificat de funcțiile Gamma, Beta și Digamma, cu idei timpurii a 20 -a din geometrie și topologie (simplex) și statistici moderne bayesiene.
(2)-simplex
Un simplex este o generalizare a noțiunii unui triunghi la mai multe dimensiuni. În mod informal, în dimensiuni k, un simplex este cel mai simplu poligon care este coca convexă a vertexurilor sale K. Vectorii care cuprind simplex trebuie să fie non-negativi și să se oprească la 1. Deci, un simplex este un mod natural de a reprezenta probabilitățile care se ridică la 1 în spații multidimensionale.
Suportul pentru distribuția Dirichlet tridimensională, punctele pe care este definită distribuția, este un (2)-simplex subsetul triunghiular al unui plan 2-dimensional care intersectează axele euclidiene la punctele (1,0,0), (0,1,0) și (0,0,1). (Orientați triunghiul în graficul interactiv de mai jos, astfel încât planul de referință să fie deasupra și vârful îndreaptă în jos și veți vedea cum se aliniază axele.)
R pachete utilizate în acest post
library(ggplot2) library(gganimate) library(dplyr) library(magick) library(MCMCpack) # for rdirichlet library(gtools) # for ddirichlet #library(patchwork) # for combining plots library(threejs) library(extraDistr)
Afișați codul
set.seed(42)
# Sample from Dirichlet distribution over 3 categories
n_samples <- 2000
alpha <- c(1, 1, 1) # uniform prior over the simplex
samples <- rdirichlet(n_samples, alpha)
# 3D coordinates: each row is (x, y, z)
x <- samples(,1)
y <- samples(,2)
z <- samples(,3)
# Visualize using threejs scatterplot
scatterplot3js(x = x, y = y, z = z,
color = "steelblue",
size = 0.2,
bg = "black",
main = "2-simplex",
axisLabels = c(
"(1,0,0)",
"(0,1,0)",
"(0,0,1)"
))
Când 
densitatea dirichlet este simetrică în mijlocul simplexului, 
. În cazul special când 
densitatea este uniformă peste simplex. Când toate 
densitatea este concentrată la vârfurile simplexului și când 
densitatea este concentrată în centrul simplexului cu cea mai mare parte a masei concentrate pe câteva valori.
Modul în care se schimbă distribuția Symetric Dirichlet ca 
modificări
Următoarea animație, care proiectează complotul de mai sus pe două dimensiuni, arată modul în care distribuția Dirichlet se schimbă ca valoare comună a 
Cod pentru funcții de ajutor
# This function generates the data frame for the animation
generate_dirichlet_animation_data <- function(alpha1_values, alpha2_values, alpha3_values, n_samples = 2000) {
library(MCMCpack)
# Triangle vertices
v1 <- c(1, 0)
v2 <- c(0, 1)
v3 <- c(0, 0)
# Projection function
project_to_triangle <- function(x1, x2, x3) {
x <- x1 * v1(1) + x2 * v2(1) + x3 * v3(1)
y <- x1 * v1(2) + x2 * v2(2) + x3 * v3(2)
data.frame(x = x, y = y)
}
# Generate animation data
animation_data <- do.call(rbind, lapply(seq_along(alpha1_values), function(i) {
alpha <- c(alpha1_values(i), alpha2_values(i), alpha3_values(i))
samples <- rdirichlet(n_samples, alpha)
projected <- project_to_triangle(samples(,1), samples(,2), samples(,3))
projected$alpha1 <- alpha(1)
projected$alpha2 <- alpha(2)
projected$alpha3 <- alpha(3)
projected$frame <- i
projected
}))
return(animation_data)
}
# This function creates the animated plot
plot_dirichlet_evolution <- function(animation_data) {
library(ggplot2)
library(gganimate)
library(grid) # for arrow units
# Triangle vertices
v1 <- c(1, 0)
v2 <- c(0, 1)
v3 <- c(0, 0)
# Label positions
label_df <- data.frame(
x = c(v1(1), v2(1), v3(1)),
y = c(v1(2), v2(2), v3(2)),
label = c("(1,0,0)", "(0,1,0)", "(0,0,1)"),
nudge_x = c(-0.08, 0.1, 0),
nudge_y = c(0, 0, 0.08)
)
# Triangle outline
triangle_df <- data.frame(
x = c(v1(1), v2(1), v3(1), v1(1)),
y = c(v1(2), v2(2), v3(2), v1(2))
)
# Build plot
p_animation <- ggplot(animation_data, aes(x = x, y = y)) +
geom_point(alpha = 0.3, size = 0.8, color = "steelblue") +
geom_density_2d(color = "red", alpha = 0.7) +
geom_path(data = triangle_df, aes(x = x, y = y), color = "black", size = 1) +
annotate("segment", x = 1, y = 0, xend = 0, yend = 1,
arrow = arrow(length = unit(0.2, "cm")),
color = "darkgreen", size = 1.2) +
xlim(0, 1) + ylim(0, 1) +
geom_text(data = label_df,
aes(x = x, y = y, label = label),
nudge_x = label_df$nudge_x,
nudge_y = label_df$nudge_y,
size = 4, fontface = "bold", color = "black") +
labs(
title = "Dirichlet Prior Evolution",
subtitle = "",
x = "Projected X",
y = "Projected Y",
caption = "Transition from uniform to concentrated distribution"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, hjust = 0.5),
plot.subtitle = element_text(size = 14, hjust = 0.5),
axis.title = element_text(size = 12),
plot.caption = element_text(size = 10, hjust = 0.5)
) +
transition_states(frame,
transition_length = 1,
state_length = 2) +
ease_aes('sine-in-out')
return(p_animation)
}
#| message: false
#| warning: false
#| code-fold: true
#| code-summary: "Animation Code"
#|
set.seed(42)
# Parameters
alpha1_values <- seq(1, 0.1, by = -0.01) # starts high, ends low
alpha2_values <- seq(1, 0.1, by = -0.01) # starts low, ends high
alpha3_values <- seq(1, 0.1, by = -0.01) # starts low, ends high
animation_data <- generate_dirichlet_animation_data(alpha1 = alpha1_values,
alpha2 = alpha2_values,
alpha3 = alpha3_values,
n_samples = 2000)
plot_dirichlet_evolution(animation_data)
Următoarele două parcele, primele și ultimele cadre ale animației, arată clar cum densitatea trece de la distribuirea uniformă a simplexului la a fi concentrat la vârfurile simplexului. La modelarea dezvoltării unui lanț Markov multi-state, așa cum făceam în postarea la care am făcut aluzie mai sus, este o practică obișnuită să selectez un Dirichlet uniform înainte de 
. Cu toate acestea, dacă credeți că procesul este probabil să înceapă distribuit uniform între state, atunci un anterior cu .png)
ar putea fi adecvat. Dacă aveți motive să credeți că procesul va începe concentrat pe anumite stări, atunci puteți explora folosind o distribuție asimetrică prin setarea valorilor diferite pentru
Afișați codul
ggplot(subset(animation_data, frame == 1), aes(x = x, y = y)) +
geom_point(alpha = 0.3, color = "darkblue") +
ggtitle("Initial Frame: Uniform Prior: alpha = (1,1,1)")
Afișați codul
ggplot(subset(animation_data, frame == max(animation_data$frame)), aes(x = x, y = y)) +
geom_point(alpha = 0.3, color = "red") +
ggtitle("Final Frame: Concetrated Prior: alpha = (.1,.1,.1)")
Această concentrare de densitate ca 
creșterea este foarte evidentă în această următoare simulare ca trece de la (0,1, 0,1, 0,1) la (10,0, 10,0, 10,0). Aici vedem distribuția concentrându -se pe medie, = (1/3, 1/3, 1/3).
Afișați codul
#Define alpha trajectories alpha1 <- seq(.1, 10, by = 0.1) alpha2 <- seq(.1, 10, by = 0.1) alpha3 <- seq(.1, 10, by = 0.1) animation_data <- generate_dirichlet_animation_data(alpha1, alpha2, alpha3, n_samples = 2000) anim2 <- plot_dirichlet_evolution(animation_data) anim2
Rețineți că animația trece prin (0,5, 0,5, 0,5), care este Jeffreys anterior pentru distribuția Dirichlet.
Varianță și entropie diferențială
Articolul Wikipedia pentru distribuția Dirichlet afișează în mod evident entropia diferențială a distribuției:
unde şi 
sunt definite mai sus și 
este funcția digamma. (Această ecuație mi -a declanșat mențiunea despre funcția digamma de mai sus.) Dar vă rugăm să vă sfătuiți ca entropia diferențială să fie definită ca: 
pentru distribuții continue nu este același lucru cu entropia Shannon 
pentru distribuții discrete și nu are o interpretare similară. Printre altele, entropia diferențială poate fi negativă, nu este invariabilă în cadrul unei schimbări de variabile și, probabil, nu se conformează nici unei intuiții pe care le -ați dezvoltat despre entropia maximă. Graficul din stânga de mai jos arată comportamentul entropiei diferențiale pentru distribuția simetrică Dirichlet pe care am avut în vedere trece de la (0,1, 0,1, 0,1) la (5,0, 5,0, 5,0). Rețineți că entropia continuă să crească dincolo de acest punct 
care corespunde distribuției uniforme pe simplex.
Afișați codul
# Define the number of dimensions for the symmetric Dirichlet distribution
n <- 3
# Define a function to calculate the entropy of a symmetric Dirichlet distribution
dirichlet_entropy <- function(alpha, n) {
# This formula is for a symmetric Dirichlet distribution
# where all parameters a_i are equal to alpha.
# The formula uses the digamma function (psi) and the log-gamma function (lgamma).
# H(alpha) = log(Beta(alpha)) + (n * alpha - n) * psi(n * alpha) - n * (alpha - 1) * psi(alpha)
# log(Beta(alpha)) = n * lgamma(alpha) - lgamma(n * alpha)
# The final formula is simplified for symmetric case.
# lgamma(x) is the log-gamma function
# digamma(x) is the psi function, the first derivative of lgamma(x)
entropy_val <- lgamma(n * alpha) - n * lgamma(alpha) +
(n * alpha - n) * digamma(n * alpha) -
n * (alpha - 1) * digamma(alpha)
return(entropy_val)
}
# Create a sequence of alpha values from 0.1 to 5
alpha_values <- seq(0.1, 5, by = 0.01)
# Calculate the entropy for each alpha value
entropy_values <- sapply(alpha_values, dirichlet_entropy, n = n)
# Create a data frame for plotting
plot_data <- data.frame(alpha = alpha_values, entropy = entropy_values)
# Plot the entropy values
ggplot(plot_data, aes(x = alpha, y = entropy)) +
geom_line(size = 1.2, color = "steelblue") +
labs(
title = "Entropy of a Symmetric Dirichlet Distribution",
subtitle = paste("for n =", n, "dimensions"),
x = "Alpha (Concentration Parameter)",
y = "Entropy"
) +
geom_vline(xintercept = 1, linetype = "dashed", color = "red", size = 1) +
annotate("text", x = 1.1, y = max(entropy_values) * 0.9,
label = "Alpha = 1 (Uniform Distribution)",
color = "red", hjust = 0, size = 4) +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
Afișați codul
k <- seq(.1, 10, by = 0.1)
mean_dir <- numeric(length(k))
var_dir <- numeric(length(k))
for (i in seq_along(k)) {
alpha <- k(i)
sum_alpha <- alpha * 3
mean_dir(i) <- alpha / sum_alpha
var_dir(i) <- (alpha * (sum_alpha - alpha)) / (sum_alpha^2 * (sum_alpha + 1))
}
df <- data.frame(k = k, mean = mean_dir, variance = var_dir)
ggplot(df, aes(x = k, y = variance)) +
geom_line(color = "blue", size = 1) +
labs(
title = "Variance of Dirichlet Distribution vs Alpha Parameter",
x = "Alpha Parameter (k)",
y = "Variance"
)
Graficul din dreapta arată variația distribuției Dirichlet ca funcție în creștere . Vedem că variația scade spre zero ca crește. Acest lucru se reflectă în a doua animație de mai sus, care arată distribuția concentrându -se pe medie. Incertitudinea va merge la zero, dar entropia diferențială se trage spre infinit. Dacă sunteți totuși intrigat de entropia diferențială, poate doriți să aruncați o privire asupra referințelor pe care le -am inclus mai jos.
