Din nou procese gaussiene | R-bloggeri

Iată problema cu care te confrunți: ai niște date zgomotoase generate de o funcție necunoscută (continuă, cu valoare reală, sexy). Nu poți accepta misterul și trebuie să cartografiezi regiunile neobservate ale funcției. În special, doriți să găsiți valoarea maximă a acestuia. Dar nu este simplu, pentru că din orice motiv funcția mister este costisitoare de rulat.

mystery_function <- function(x) {
  cost(dollars = 100)
  1.1 * sin(sqrt(x))
}

Imaginați-vă că sunt taxele API, consumul de simboluri sau tarife sau ceva de genul. Ideea este că nu poți să-ți faci drum prin forța brută prin asta. Ai nevoie de ajutor.

Iată observațiile pe care le aveți până acum. E puțin zgomot și pe ei, doar pentru a îngreuna viața.

Cod

domain <- c(0, 40)

observations <- tibble(x = c(0, 8, 15, 20, 35)) |>
  mutate(
    y = mystery_function(x) + rnorm(length(x), mean = 0, sd = 0.1),
  )

ggplot(observations) +
  aes(x = x, y = y) +
  scale_x_continuous(limits = domain) +
  geom_point(colour = "red", size = 3) +
  geom_function(fun = mystery_function, colour = "black", linetype = "dashed")

Ai putea să te potrivească doar cu o linie, sigur. Dacă aveți o idee bună despre funcția care vă produce datele, acesta este cel mai simplu și mai evident lucru de făcut.

Cod

ggplot(observations) +
  aes(x = x, y = y) +
  scale_x_continuous(limits = domain) +
  geom_point(colour = "red", size = 3) +
  geom_function(fun = mystery_function, colour = "black", linetype = "dashed") +
  geom_smooth(method = "lm", formula = y ~ sin(x**0.5))

Dacă nu am fi stăpâni în acest univers de jucării, s-ar putea să nu avem idee ce funcție să ne potrivim. Am putea încadra orice altceva, ca un polinom. Dar orice ne potrivim, facem o mulțime de presupuneri despre funcție și ar putea ajunge destul de departe de marcaj, mai ales dacă nu avem multe observații.

Cod

ggplot(observations) +
  aes(x = x, y = y) +
  scale_x_continuous(limits = domain) +
  geom_point(colour = "red", size = 3) +
  geom_function(fun = mystery_function, colour = "black", linetype = "dashed") +
  geom_smooth(method = "lm", formula = "y ~ poly(x, 3)")

Ceea ce am făcut până acum este abordarea non-bayesiană a problemei. De ce să vă deranjați deloc cu abordarea bayesiană (procesele gaussiene)?

• În primul rând, pentru că nu avem neapărat o idee bună care este funcția noastră. Ar putea fi foarte neliniar. Un proces gaussian este un model neparametric, adică nu presupune nicio formă specifică a funcției.

• În al doilea rând, ca de obicei cu modelele bayesiene, puteți obține o înțelegere mai bogată a incertitudinii pe care o puteți obține cu modele liniare și intervale de încredere. De asemenea, este mai potrivit pentru scenariile cu date reduse.

Totuși, să-ți înțelegi abordarea bayesiană poate fi puțin dificil. Trebuie să uităm de ideea de a optimiza coeficienții unei formule liniare pentru a se potrivi modelului nostru. În schimb, estimăm posibilele locații ale liniei. Nu ne mai pasă de formulă.

Mai degrabă decât o formulă, putem folosi un proces Gaussian (GP). GP este o distribuție care generează „funcții” (linii). Cum așa? Este o distribuție normală multivariată, care poate genera o valoare aleatorie pentru mai multe dimensiuni.

Dacă vă imaginați că fiecare dimensiune este un punct pe axa x, un eșantion din distribuția normală multivariată descrie o linie aleatorie sau o funcție necunoscută, pentru a o privi în alt mod.

Cod

n.dims <- 16

# Sigma specifies the covariance matrix of the variables.
# This identify matrix means each dim is completely unrelated.
sigma <- matrix(0, nrow = n.dims, ncol = n.dims)
for (i in 1:n.dims) {
  sigma(i, i) = 1
}

random.draw <- MASS::mvrnorm(
  n = 1, # n.samples
  mu = rep(0, n.dims), # all centered on 0
  Sigma = sigma
)

ggplot(
  data.frame(x = 1:n.dims, y = random.draw)
) +
  aes(x = x, y = y) +
  geom_line()

The nucleu (sigma în codul de mai sus) controlează forma funcțiilor generate. În exemplul de mai sus este o matrice de identificare, dar o alegere mai convențională ar fi covarianța exponențială la pătrat. Acest lucru generează o varietate frumoasă de funcții. Cei doi hiperparametri alfa și lambda controlează ceva precum amplitudinea și frecvența funcțiilor generate.

Cod

# Kernel
sq_exp_cov <- function(x, lambda, alpha) {
  n <- length(x)
  K <- matrix(0, n, n)
  for (i in 1:n) {
    for (j in 1:n) {
      diff <- sqrt(sum((x(i) - x(j))^2))
      K(i, j) <- alpha^2 * exp(-diff^2 / (2 * lambda^2))
    }
  }
  K
}

plot_random_gp <- function(lambda, alpha) {
  n.samples <- 5
  x <- seq(0, 9, 0.1)
  samples <- MASS::mvrnorm(
    n = n.samples,
    mu = rep(0, length(x)),
    Sigma = sq_exp_cov(x, lambda = lambda, alpha = alpha)
  )

  t(samples) |>
    as_tibble() |>
    mutate(x = x) |>
    pivot_longer(-x, names_to = "sample", values_to = "y") |>
    ggplot() +
    aes(x = x, y = y, group = sample, colour = sample) +
    geom_line()
}
plot_random_gp(lambda = 2, alpha = 1)

Deci avem acest generator de funcții misterioase! Poate propune funcții posibile și putem vedea cât de bine se potrivesc cu datele. De asemenea, putem prezice valori pentru aceste funcții în puncte pentru care nu avem observații.

E timpul pentru Stan. Conform ultimului blog, acest cod Stan a fost preluat dintr-un curs excelent de Tâmplărie.

Susțin că cel mai bun mod de a-l citi pe Stan este înapoi în față, așa că iată câteva Stan invalide pe care l-am inversat și adnotat.

```{stan}
model {
  // f_x is sampled from the GP.
  // Remember it's a vector of f(x), not parameters describing f.
  f_x ~ multi_normal(rep_vector(0, n), K);

  // Likelihood is evaluated on f(x) plus normally distributed noise.
  y_obs ~ normal(f(1 : n_obs), sigma);

}
parameters {
  // Our outputs will be f(x) for all x (our observations, x_obs, and the predictions we want, x_pred)
  vector(n) f_x;
}
data {
  // Observed data
  int n_obs;
  array(n_obs) real x_obs;
  array(n_obs) real y_obs;
  
  // Observation error - pick a value
  real sigma;
  
  // x values for which we aim to predict f(x)
  int n_pred;
  array(n_pred) real x_pred;
  
  // Hyperparameters for the kernel
  real alpha;
  real lambda;
}
transformed data {
  // We join the x observations and desired x prediction points
  // because we want f(x) for both observed data and new predictions
  int n = n_obs + n_pred;
  
  array(n) real x;
  x(1 : n_obs) = x_obs;
  x((n_obs + 1): n) = x_pred;
  
  // We calculate the Kernel values for all observed x
  matrix(n, n) K;
  K = gp_exp_quad_cov(x, alpha, lambda);
  
  // Add "nugget" on diagonal for numerical stability
  for (i in 1 : n) {
    K(i, i) = K(i, i) + 1e-6;
  }
}

```

Fișierul Stan valid este aici. O versiune optimizată cu descompunere Cholesky este aici, care este cea pe care o vom rula.

x.pred <- seq(domain(1), domain(2), by = 0.25)
x.vals <- c(observations$x, x.pred)

model <- cmdstan_model(stan_file = "gp-cholesky.stan", exe = "gp.stan.bin")
sample <- model$sample(
  seed = SEED,
  list(
    n_data = nrow(observations),
    x_data = as.array(observations$x),
    y_data = as.array(observations$y),
    # Sigma here is our guess about the noise level.
    sigma = 0.2,
    n_pred = length(x.pred),
    x_pred = x.pred,
    alpha = 1,
    lambda = 2
  ),
  parallel_chains = 4,
  max_treedepth = 20,
  show_messages = FALSE # disabled to avoid polluting the blog post, should be TRUE
)

sample

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
   lp__   -83.63 -83.23 8.92 8.77 -99.08 -69.46 1.00     1503     2005
   eta(1)   0.13   0.13 0.19 0.18  -0.19   0.45 1.00     6394     2417
   eta(2)   0.27   0.27 0.19 0.20  -0.05   0.58 1.00     6378     2810
   eta(3)  -0.67  -0.68 0.20 0.20  -1.00  -0.35 1.00     6100     3056
   eta(4)  -0.94  -0.94 0.20 0.21  -1.27  -0.62 1.00     5973     2324
   eta(5)  -0.34  -0.34 0.20 0.19  -0.66  -0.02 1.00     6341     2760
   eta(6)  -0.01  -0.02 1.03 1.02  -1.69   1.69 1.00     6262     2636
   eta(7)  -0.01  -0.02 0.99 0.99  -1.65   1.65 1.00     6661     2761
   eta(8)   0.01   0.03 0.97 0.99  -1.59   1.62 1.00     5957     3165
   eta(9)   0.00   0.01 1.02 1.02  -1.69   1.66 1.00     5475     2417

 # showing 10 of 333 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

Rhat și ESS sunt sănătoase (≤1,01 și respectiv >>100), așa că putem fi fericiți. Regrupăm eșantioanele în valori (x, y) și le trasăm alături de observațiile originale. Acum arată un pic ca un decor de Crăciun.

Cod

tidy_sample <- function(sample) {
  mat <- sample$draws(format = "draws_matrix")
  as_tibble(mat) |>
    # Every row is a draw, which we number
    mutate(draw = 1:n()) |>
    # Since each column is an observation f(x) for x indices, we pivot
    pivot_longer(starts_with("f"), names_to = "x", values_to = "y") |>
    # And map the x index back to an x value
    mutate(
      idx = as.numeric(str_extract(x, "(0-9)+")),
      x = x.vals(idx),
      y = as.numeric(y)
    )
}

plot_draws <- function(sample) {
  draws <- tidy_sample(sample)
  ggplot(observations) +
    aes(x = x, y = y) +
    geom_line(
      data = draws,
      mapping = aes(group = draw),
      alpha = 0.01,
      colour = "darkgreen"
    ) +
    geom_point(colour = "red", size = 3) +
    geom_function(fun = mystery_function, colour = "black", linetype = "dashed")
}

plot_draws(sample)

Modelul bayesian are un model de incertitudine mult mai bogat decât potrivirea liniară. De exemplu, este mult mai puțin incert în ceea ce privește pozițiile în care avem date.

Vă amintiți obiectivul nostru inițial de a găsi valoarea maximă a funcției mister? L-am putea restrânge cu mai multe date. Cel mai util ar fi dacă am avea o modalitate de a determina următoarea cea mai bună observație de obținut. Să presupunem că căutăm valoarea maximă a . Dacă am mai adăuga o observație la valoarea alegerii noastre, ce ar trebui asta fi?

Acesta este functia de achizitiepe care îl maximizăm pentru a găsi cel mai bun în continuare a dobândi. Există două obiective de echilibrat: unde este posibilă îmbunătățirea și unde există o îmbunătățire potențială mare.

Cea mai populară funcție de achiziție este Expected Improvement, care face exact asta. Pentru fiecare calculează valoarea așteptată a îmbunătățirii.

Deoarece avem un model bayesian cu multe mostre posterioare, matematica este destul de ușoară. Ne uităm doar la fiecare punct și luați îmbunătățirea medie față de cele mai bune (adică maxim dacă ne propunem să maximizăm ). Oricare ar fi a dat cea mai mare îmbunătățire medie este norocosul câștigător și ar trebui să fie obținut în continuare în experiment.

y.best <- max(observations$y)

expected_improvement <- function(sample, x.vals, y.best) {
  mat <- sample$draws(format = "draws_matrix")
  mat <- mat(, str_starts(colnames(mat), "f"))
  mat <- t(mat)

  n.points <- nrow(mat)
  ei <- numeric(length = n.points)

  for (i in 1:n.points) {
    samples <- mat(i, )
    improvements <- pmax(0, samples - y.best)
    ei(i) <- mean(improvements)
  }

  x.vals(which.max(ei))
}

x.next <- expected_improvement(sample, x.vals, y.best)
sprintf("Next x: %.2f", x.next)

Și voilà, modelul spune că ar trebui să măsurăm la x=5. Maximul funcției mister atunci când este rezolvată analitic este de fapt de aproximativ 2,47, deci acesta este intuitiv un pas în direcția corectă. Nu-i rău, funcție de achiziție!

Dacă măsurăm din nou, ce atunci? Am putea clăti și repeta, până ajungem la maxim. În esență, asta este tot ceea ce înseamnă optimizarea bayesiană: potriviți modelul, calculați următoarea achiziție cu funcția de îmbunătățire așteptată și potriviți din nou.

Desigur, comunitatea R dă; puteți folosi tune_bayes în loc să rulați propriul optimizator.

Deci, iată-l. Am adaptat un model Gaussian Process pe observațiile zgomotoase dintr-o funcție misterioasă și am folosit funcția de achiziție a îmbunătățirii așteptate pentru a determina următoarea observație pe care ar trebui să o obținem dacă dorim să găsim maximul funcției, care este primul pas în optimizarea bayesiană.

Aștept cu nerăbdare să te revăd peste șase luni când am uitat toate astea și trebuie să le reînvăț încă o dată!