După ce a trecut prin evaluarea inter pares rOpenSci, pachetul dicționar este acum disponibil nou pe CRAN. Vă permite să faceți rapid distribuții de probabilitate – fie din câteva intrări, fie din biblioteca sa încorporată – și apoi să le examinați în detaliu.
Aceste distribuții formează blocurile de bază care împletesc modelele statistice avansate cu ecosistemul probaver mai larg, care este construit pentru a elibera modelatorii de codificarea la nivel scăzut, astfel încât conductele de producție să rămână prietenoase cu oamenii. Momentan, celelalte pachete probaverse sunt distplyr, permițându-vă să transformați distribuțiile în forme noi și famish, permițându-vă să reglați distribuțiile la date. Dezvoltate ținând cont de cazuri de utilizare a analizei de risc, cum ar fi clima și asigurările, aceleași instrumente se traduc fără probleme în simulări, predare și alte setări aplicate.
Această postare evidențiază primele 3 caracteristici ale acestei cele mai tinere versiuni de dicționar. Să începem prin a încărca pachetul.
library(distionary)
Caracteristica 1: mai mult decât distribuții Base R
Desigur, toate distribuțiile Base R sunt disponibile în dicționar. Iată distribuția normală preferată a tuturor.
dst_norm(0, 1) Normal distribution (continuous) --Parameters-- mean sd 0 1 plot(dst_norm(0, 1))
Și bunul Poisson.
dst_pois(3) Poisson distribution (discrete) --Parameters-- lambda 3 plot(dst_pois(3))
Dar sunt incluse și distribuții suplimentare care schimbă jocul.
O Distribuție nulăcare evaluează întotdeauna la NA. Când rulați un algoritm care întâmpină o problemă, puteți returna o distribuție Nulă în loc să aruncați o eroare. Nici măcar pașii de evaluare din aval nu vor eroa, deoarece codul încă vede o distribuție mai degrabă decât un gol NA sau NULL.
# Make a Null distribution. null <- dst_null() # Null distributions always evaluate to NA. eval_quantile(null, at = c(0.25, 0.5, 0.75)) (1) NA NA NA mean(null) (1) NA
Distribuții empiriceunde datele sunt distributia. Acestea respectă comportamentul observat fără a forța o formă specifică și sunt, de asemenea, utilizate în mod obișnuit ca reper pentru comparație cu alte modele. Iată un exemplu folosind concentrația de ozon din airquality setul de date care vine încărcat cu R.
# Empirical distribution of Ozone from the `airquality` dataset. emp <- dst_empirical(airquality$Ozone, na_action_y = "drop") # Inspect print(emp, n = 5) Finite distribution (discrete) --Parameters-- # A tibble: 67 × 2 outcomes probs1 1 0.00862 2 4 0.00862 3 6 0.00862 4 7 0.0259 5 8 0.00862 # ℹ 62 more rows
Comparați funcția sa de distribuție cumulativă (CDF) cu cea a unei distribuții Gamma adaptate nivelurilor de ozon, împrumutând pachetul de foame al probaversei pentru sarcina de adaptare.
# Fit a Gamma distribution to Ozone using the famish package.
library(famish)
gamma <- fit_dst_gamma(airquality$Ozone, na_action = "drop")
# Plot the cumulative distribution functions (CDFs) together.
plot(emp, "cdf", n = 1000, xlab = "Ozone Levels (ppb)")
plot(gamma, "cdf", add = TRUE, col = "red")
legend(
"bottomright",
legend = c("Empirical", "Fitted Gamma"),
col = c("black", "red"),
lty = 1
)
Aceste distribuții de manuale devin mult mai utile odată ce devin blocuri pentru construirea unui sistem. De exemplu, ar putea forma distribuții predictive într-un context de învățare automată sau pot fi legate de alte variabile. Aceasta este ceea ce probaverse caută să facă posibil.
Caracteristica 2: prietenos cu fluxurile de lucru tabulare ordonate
Mai întâi, încărcați tidyverse pentru a activa fluxuri de lucru tabelare ordonate. Și da, probaverse este numit după tidyverse deoarece își propune să fie un „tidyverse pentru probabilitate”.
library(tidyverse)
Puteți ignora cu siguranță această bucată următoare, dacă nu doriți să vedeți cum mă disput niște date financiare pentru dvs.
# Wrangle the stocks data frame using tidyverse. stocks <- as_tibble(EuStockMarkets) |> mutate(across(everything(), (x) 100 * (1 - x / lag(x)))) |> drop_na()
The stocks datele pe care le-am disputat sunt un tabel pierdere procentuală zilnică pentru patru mari indici bursieri europeni. Datele nu contează pentru acest exemplu, așa că au fost omise.
stocks # A tibble: 1,859 × 4 DAX SMI CAC FTSE1 0.928 -0.620 1.26 -0.679 2 0.441 0.586 1.86 0.488 3 -0.904 -0.328 0.576 -0.907 4 0.178 -0.148 -0.878 -0.579 5 0.467 0.889 0.511 0.720 6 -1.25 -0.676 -1.18 -0.855 7 -0.578 -1.23 -1.32 -0.824 8 0.287 0.358 0.193 -0.0837 9 -0.637 -1.11 -0.0171 0.522 10 -0.118 -0.437 -0.314 -1.41 # ℹ 1,849 more rows
În primul rând, să ne concentrăm asupra indicelui bursier DAX. Potriviți o distribuție empirică ca data trecută (observați că folosesc o mască de date în dst_empirical() de data asta).
# Fit an empirical distribution to the DAX stock index. dax <- dst_empirical(DAX, data = stocks) # Inspect the CDF. plot(dax, xlab = "Daily Loss (%)")
Puteți calcula cu ușurință unele cuantile standard în format tabelar, astfel încât intrările să fie plasate alături de ieșirile calculate: utilizați doar enframe_ prefix în loc de eval_ așa cum am făcut mai sus cu distribuția Null.
enframe_quantile(dax, at = c(0.25, 0.5, 0.75), arg_name = "prob") # A tibble: 3 × 2 prob quantile1 0.25 -0.638 2 0.5 -0.0473 3 0.75 0.468
Sau, mai la obiect aici – și atrăgând la punctul slab al probaverse pentru munca axată pe risc – puteți calcula niveluri de rentabilitate (cunoscute și sub numele de „Valoare la risc” în aplicațiile financiare) pentru anumite perioade de rentabilitate. Dacă nu știți care sunt acestea, sunt doar nume de lux pentru cuantile.
return_periods <- c(5, 50, 100, 200, 500) enframe_return( dax, at = return_periods, arg_name = "return_period", fn_prefix = "daily_loss_pct" ) # A tibble: 5 × 2 return_period daily_loss_pct1 5 0.621 2 50 2.17 3 100 2.75 4 200 3.08 5 500 3.71
Ieșirea tabelară devine și mai puternică atunci când este inserată într-un tabel de modele, deoarece facilitează comparațiile și tendințele. Pentru a demonstra, construiți un model pentru fiecare stoc. Mai întâi, prelungiți datele pentru această sarcină.
# Lengthen the data using tidyverse. stocks2 <- pivot_longer( stocks, everything(), names_to = "stock", values_to = "daily_loss_pct" ) # Inspect stocks2 # A tibble: 7,436 × 2 stock daily_loss_pct1 DAX 0.928 2 SMI -0.620 3 CAC 1.26 4 FTSE -0.679 5 DAX 0.441 6 SMI 0.586 7 CAC 1.86 8 FTSE 0.488 9 DAX -0.904 10 SMI -0.328 # ℹ 7,426 more rows
Construiți un model pentru fiecare stoc folosind a group_by + summarise flux de lucru din tidyverse (vă rugăm să scuzați nevoia actuală de a încheia distribuția în list()). Observați că distribuțiile devin intrări de tabel, indicate aici de clasa lor
# Create an Empirical distribution for each stock. models <- stocks2 |> group_by(stock) |> summarise(model = list(dst_empirical(daily_loss_pct))) # Inspect models # A tibble: 4 × 2 stock model1 CAC
2 DAX 3 FTSE 4 SMI
Acum puteți utiliza un flux de lucru ordonat pentru a calcula tabele de cuantile pentru fiecare model și pentru a le extinde. De fapt, acest flux de lucru este suficient de comun încât mă gândesc să adaug un verb dedicat pentru el.
return_levels <- models |> mutate( df = map( model, enframe_return, at = return_periods, arg_name = "return_period", fn_prefix = "daily_loss_pct" ) ) |> unnest(df) |> select(!model) # Inspect print(return_levels, n = Inf) # A tibble: 20 × 3 stock return_period daily_loss_pct1 CAC 5 0.757 2 CAC 50 2.37 3 CAC 100 2.78 4 CAC 200 3.41 5 CAC 500 3.97 6 DAX 5 0.621 7 DAX 50 2.17 8 DAX 100 2.75 9 DAX 200 3.08 10 DAX 500 3.71 11 FTSE 5 0.542 12 FTSE 50 1.58 13 FTSE 100 2.05 14 FTSE 200 2.31 15 FTSE 500 2.87 16 SMI 5 0.552 17 SMI 50 2.03 18 SMI 100 2.52 19 SMI 200 2.91 20 SMI 500 3.55
Rezultatul este un set de date ordonat, care este gata pentru majoritatea analizelor. De exemplu, puteți realiza cu ușurință o comparație a nivelurilor de rentabilitate ale fiecărui stoc. Fac aceste comploturi tot timpul pentru a facilita luarea deciziilor bazate pe riscuri.
return_levels |>
mutate(stock = fct_reorder2(stock, return_period, daily_loss_pct)) |>
ggplot(aes(return_period, daily_loss_pct, colour = stock)) +
geom_point() +
geom_line() +
theme_bw() +
scale_x_log10(
"Return Period (days)",
minor_breaks = c(1:10, 1:10 * 10, 1:10 * 100)
) +
scale_y_continuous("Daily Loss", label = scales::label_number(suffix = "%")) +
annotation_logticks(side = "b") +
scale_colour_discrete("Stock Index")
Caracteristica 3: faceți distribuția de care aveți nevoie
Puteți crea propriile distribuții cu dicționar specificând doar un set minim de proprietăți; toate celelalte proprietăți sunt derivate automat și pot fi recuperate atunci când este necesar.
Să presupunem că aveți nevoie de o distribuție Gamma inversă, dar nu este disponibilă în dicționar. În prezent, dicționarul presupune că veți furniza cel puțin densitatea și CDF; le puteți prelua din pachetul extraDistr (funcții dinvgamma() şi pinvgamma()). Conectați-le la dicționare distribution() funcția și bucurați-vă de acces la o varietate de proprietăți pe care nu le-ați specificat, cum ar fi funcția medie, varianță, asimetrie și pericol.
# Make an Inverse Gamma distribution (minimal example). ig <- distribution( density = function(x) extraDistr::dinvgamma(x, alpha = 5, beta = 20), cdf = function(x) extraDistr::pinvgamma(x, alpha = 5, beta = 20), .vtype = "continuous", ) # Calculate anything. mean(ig) (1) 5 variance(ig) (1) 8.333333 skewness(ig) (1) 3.464085 plot(ig, "hazard", to = 20, n = 1000, xlab = "Outcome")
De asemenea, ați putea lua în considerare distribuirea a .name – se plătește atunci când jonglați cu mai multe modele. Adăugând .parameters oferă o specificitate suplimentară distribuției .namedar nu sunt încă utilizate în scopuri funcționale.
Iată o implementare mai completă a distribuției Inverse Gamma, de data aceasta implementată în funcție de cei doi parametri. Observați că verific, de asemenea, dacă parametrii sunt pozitivi ( salutări pentru pachetul de șahmat).
dst_invgamma <- function(alpha, beta) {
checkmate::assert_number(alpha, lower = 0)
checkmate::assert_number(beta, lower = 0)
distribution(
density = (x) extraDistr::dinvgamma(x, alpha = alpha, beta = beta),
cdf = (x) extraDistr::pinvgamma(x, alpha = alpha, beta = beta),
quantile = (p) extraDistr::qinvgamma(p, alpha = alpha, beta = beta),
random = (n) extraDistr::rinvgamma(n, alpha = alpha, beta = beta),
.name = "Inverse Gamma",
.parameters = list(alpha = alpha, beta = beta),
.vtype = "continuous",
)
}
Acum putem face aceeași distribuție Gamma inversă ca înainte:
ig2 <- dst_invgamma(5, 20) # Inspect ig2 Inverse Gamma distribution (continuous) --Parameters-- alpha beta 5 20
Apropo, această caracteristică – posibilitatea de a inspecta alte proprietăți de distribuție chiar și atunci când acestea nu sunt specificate – este excelentă pentru a afla despre probabilitate. Acest lucru se datorează faptului că puteți vedea numeroasele moduri în care distribuțiile pot fi reprezentate, nu doar prin funcțiile obișnuite de densitate sau de masă de probabilitate văzute în manuale.
Această caracteristică permite, de asemenea, extensibilitatea probaversei. De exemplu, pachetul distplyr al probaversei creează distribuții de amestec, care nu au o formulă explicită pentru funcția cuantilă. Cu toate acestea, acest lucru nu este problematic – distribuția poate fi încă definită, iar dicționarul va da seama care sunt cuantilele.
Ce urmează?
În prezent, pachetul dicționar oferă funcționalitate cheie pentru definirea și evaluarea obiectelor de distribuție. Obiectivele viitoare includ:
- Acoperire mai largă. Extinderea dincolo de distribuțiile continue univariate, astfel încât cazurile de utilizare mixte discrete/continue și multivariate se simt de primă clasă.
- Valori native de risc. Efectuarea funcțiilor de cost, a așteptărilor finale și a altor valori de decizie să fie evaluate mai natural.
- Ergonomia fluxului de lucru. Relaxarea cerinței de a furniza perechi densitate/CDF și adăugarea de verbe care simplifică comunitatea
mutate–map–unnestmodele.
Dacă acest lucru te entuziasmează, alătură-te conversației deschizând o problemă sau contribuind.
Mulțumiri speciale recenzenților rOpenSci Katrina Brock și Christophe Dutang pentru comentariile perspicace care au îmbunătățit acest pachet. De asemenea, mulțumesc lui BGC Engineering Inc., Consorțiului R și Agenției Spațiale Europene împreună cu Politecnico di Milano pentru sprijinirea acestui proiect.
