Saïd Maanan's questions

Tenho a seguinte carga elétrica diária para a Espanha:

> head(spain_daily)
# A tibble: 6 × 2
  Date       Daily_Load
  <date>          <dbl>
1 2022-01-01     20872.
2 2022-01-02     22729.
3 2022-01-03     27358.
4 2022-01-04     28262.
5 2022-01-05     28048.
6 2022-01-06     24182.
> tail(spain_daily)
# A tibble: 6 × 2
  Date       Daily_Load
  <date>          <dbl>
1 2024-12-26     25210.
2 2024-12-27     26363.
3 2024-12-28     25041.
4 2024-12-29     24182.
5 2024-12-30     27073.
6 2024-12-31     25924.

Como o consumo elétrico é cíclico, decompus a série usando os seguintes comandos:

library(forecast)
library(stR)

autoplot(mstl(msts(spain_daily$Daily_Load, seasonal.periods = c(7, 365.25))))

O que me deu o seguinte enredo:

Tentei recriar o gráfico acima usando a função stRdo pacote AutoSTRcom o seguinte comando:

plot(AutoSTR(msts(spain_daily$Daily_Load, seasonal.periods = c(7, 365.25))))

E me dá o seguinte gráfico:

Como tornar este gráfico mais semelhante ao primeiro, ou seja, exibir a série de tendências em seu próprio gráfico e remover a série de previsão (azul)?

Obrigado.

Tenho uma lista de dados (veja abaixo), que represento graficamente usando este script:

# Define plotting parameters
pch_values <- c(16, 17)  # Different symbols for beta_0 and beta_1
colors <- c("red", "blue", "gray60")  # Colors for different noise types
lty_values <- c(1, 2, 3)  # Line types for different noise types
quartile_titles <- c("Absolute Bias of First Quartile", "Absolute Bias of Second Quartile (Median)", "Absolute Bias of Third Quartile")

# Directory for saving plots
output_dir <- "Abs_Quartile_Bias_Plots"
if (!dir.exists(output_dir)) dir.create(output_dir)

# Compute unified y-axis limits for all quartile plots in test_data
y_values_all <- unlist(sapply(test_data, function(noise_data) {
  sapply(noise_data, function(T_data) T_data)
}), use.names = FALSE)

ylim_min <- min(y_values_all, na.rm = TRUE)
ylim_max <- max(y_values_all, na.rm = TRUE)
margin <- 0.05 * (ylim_max - ylim_min)  # Add 5% margin above and below
unified_ylim <- c(ylim_min - margin, ylim_max + margin)

# Set output file
png(file.path(output_dir, "Abs_Quartile_Bias_Case1.png"), 
    width = 1200, height = 400, res = 120)
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))  # 3 side-by-side subplots

# Plot each absolute quartile bias
for (quartile_idx in 1:3) {
  # Start a blank plot with unified y-axis limits
  plot(NULL, log = "x", xlim = range(sample_sizes), ylim = unified_ylim, 
       xlab = "Sample Size (T)", 
       ylab = quartile_titles[quartile_idx],
       main = quartile_titles[quartile_idx])
  
  # Add lines and points for beta_0 and beta_1 for each noise type
  for (noise_idx in seq_along(noise_types)) {
    noise <- noise_types[noise_idx]
    for (coef_name_idx in seq_along(pch_values)) {
      coef_name <- c("beta_0", "beta_1")[coef_name_idx]  # Explicitly specify coefficient names
      
      y_values <- sapply(test_data[[noise]][[coef_name]], function(T_data) T_data[[quartile_idx]])
      
      # Plot lines
      lines(sample_sizes, y_values, col = colors[noise_idx], lty = lty_values[noise_idx], lwd = 1)
      # Plot points
      points(sample_sizes, y_values, col = colors[noise_idx], pch = pch_values[coef_name_idx])
    }
  }
  
  # Add legends
  if (quartile_idx == 1) {  # Only add legends to the first plot
    legend("topright", 
           legend = c(expression(beta[0]), expression(beta[1])),
           pch = pch_values, 
           col = "black", bty = "n")
    legend("topright", 
           legend = str_to_title(noise_types), 
           col = colors, 
           lty = lty_values, 
           lwd = 1, bty = "n", inset = c(.18, 0))
  }
}

dev.off()

O que produz o seguinte valor:

O problema, como você pode ver, é que ele ylimnão é ajustado para cada subtrama, o que cria gráficos em que a extremidade superior ylimfica longe do maior valor das linhas dentro do subtrama.

Pergunta:

Alguém pode me ajudar a modificar o código para que ele ylimseja ajustado dinamicamente ao intervalo de dados em questão?

Peço desculpas se meus dados não estiverem organizados de forma ideal.

Aqui está a dput()saída dos dados em questão:

> dput(test_data)
list(gaussian = list(beta_0 = list(`50` = c(`25%` = 0.243594614710333, 
`50%` = 0.00277753736411457, `75%` = 0.245669686373254), `100` = c(`25%` = 0.171918499273751, 
`50%` = 9.39373213166839e-05, `75%` = 0.183152746837558), `200` = c(`25%` = 0.116881126444175, 
`50%` = 0.00117104968639603, `75%` = 0.123610284184898), `500` = c(`25%` = 0.0804521896728776, 
`50%` = 0.00103027837221392, `75%` = 0.0784101295289328), `1000` = c(`25%` = 0.0541010700519603, 
`50%` = 0.00294703982537503, `75%` = 0.0531978110356446), `2000` = c(`25%` = 0.0377543581178926, 
`50%` = 0.00333965713516848, `75%` = 0.0413079011502395), `5000` = c(`25%` = 0.0217425538489762, 
`50%` = 0.000460109571013723, `75%` = 0.0257457441062152)), beta_1 = list(
    `50` = c(`25%` = 0.00882506794511806, `50%` = 0.000371504996548033, 
    `75%` = 0.00765082797750871), `100` = c(`25%` = 0.00290263844961158, 
    `50%` = 6.82498075157412e-05, `75%` = 0.00286586470881645
    ), `200` = c(`25%` = 0.00109095625190614, `50%` = 5.22460831402505e-05, 
    `75%` = 0.00109448558611547), `500` = c(`25%` = 0.00026571655332952, 
    `50%` = 8.85613155698906e-07, `75%` = 0.000259202833770455
    ), `1000` = c(`25%` = 9.21700578033757e-05, `50%` = 7.24859772205377e-07, 
    `75%` = 9.59841859602406e-05), `2000` = c(`25%` = 3.41775667214161e-05, 
    `50%` = 3.30163678108342e-06, `75%` = 2.86060854759462e-05
    ), `5000` = c(`25%` = 8.28928386886751e-06, `50%` = 4.18267251500737e-07, 
    `75%` = 8.43229566438453e-06))), laplace = list(beta_0 = list(
    `50` = c(`25%` = 0.198717632408904, `50%` = 0.0081310174987681, 
    `75%` = 0.210213471331161), `100` = c(`25%` = 0.130486991150159, 
    `50%` = 0.00513011870745728, `75%` = 0.151674765218005), 
    `200` = c(`25%` = 0.111212423845099, `50%` = 0.00734704246964579, 
    `75%` = 0.0981520122896047), `500` = c(`25%` = 0.0630465314163474, 
    `50%` = 0.00897075972839323, `75%` = 0.0683817884348215), 
    `1000` = c(`25%` = 0.042034483419907, `50%` = 0.00132855758101336, 
    `75%` = 0.0425837417269435), `2000` = c(`25%` = 0.0339129308142302, 
    `50%` = 0.00345531227818663, `75%` = 0.0280143719137393), 
    `5000` = c(`25%` = 0.0198791807771237, `50%` = 0.000479565883414246, 
    `75%` = 0.0194025102135911)), beta_1 = list(`50` = c(`25%` = 0.00740494352386101, 
`50%` = 0.00058286479518399, `75%` = 0.00741009508728219), `100` = c(`25%` = 0.00244184197921071, 
`50%` = 0.000134746791191187, `75%` = 0.00244511866997321), `200` = c(`25%` = 0.000866632964359182, 
`50%` = 8.7678493771115e-05, `75%` = 0.0010573916668597), `500` = c(`25%` = 0.000235006834866658, 
`50%` = 9.75544170911391e-06, `75%` = 0.000212441344117131), 
    `1000` = c(`25%` = 8.00301522252411e-05, `50%` = 1.07320472242378e-06, 
    `75%` = 7.29608951477445e-05), `2000` = c(`25%` = 2.54472054908028e-05, 
    `50%` = 1.46816056423305e-06, `75%` = 2.84447257516973e-05
    ), `5000` = c(`25%` = 6.16220803828504e-06, `50%` = 1.48241995123755e-07, 
    `75%` = 6.65929311072233e-06))), cauchy = list(beta_0 = list(
    `50` = c(`25%` = 0.274524411682779, `50%` = 0.0361625726093231, 
    `75%` = 0.336776027564404), `100` = c(`25%` = 0.198215929997078, 
    `50%` = 0.00328917426856501, `75%` = 0.207109168934052), 
    `200` = c(`25%` = 0.155318692086412, `50%` = 0.00672679937945397, 
    `75%` = 0.154703229084234), `500` = c(`25%` = 0.0880029220252172, 
    `50%` = 0.00278228951379011, `75%` = 0.0884479446139164), 
    `1000` = c(`25%` = 0.0688536372229773, `50%` = 0.00305819621170977, 
    `75%` = 0.0644673760777885), `2000` = c(`25%` = 0.0507696330400698, 
    `50%` = 0.000188404592404767, `75%` = 0.0500372607497894), 
    `5000` = c(`25%` = 0.0310586078673812, `50%` = 0.0017270922827648, 
    `75%` = 0.0322962294132709)), beta_1 = list(`50` = c(`25%` = 0.0120626033466085, 
`50%` = 0.0017021299049178, `75%` = 0.00978118814587825), `100` = c(`25%` = 0.00338466332476295, 
`50%` = 7.46446842248005e-06, `75%` = 0.0037767462745677), `200` = c(`25%` = 0.00125677406659053, 
`50%` = 4.41521154233016e-05, `75%` = 0.00133801259948108), `500` = c(`25%` = 0.000315465752873667, 
`50%` = 1.30850267132665e-05, `75%` = 0.0003114951811094), `1000` = c(`25%` = 0.000119877019459924, 
`50%` = 5.46893714759022e-07, `75%` = 0.000117584995672271), 
    `2000` = c(`25%` = 4.14058369635484e-05, `50%` = 9.97860400531181e-07, 
    `75%` = 4.36402038799244e-05), `5000` = c(`25%` = 9.47781604621056e-06, 
    `50%` = 9.22661139490799e-07, `75%` = 1.12907100797699e-05
    ))))

Estou trabalhando em um estudo de simulação em R, onde estou plotando curvas de viés e variância para diferentes estimadores. Quero que os títulos dos subplots incluam expressões matemáticas formatadas corretamente, especificamente \beta_0 e \beta_1, usando expression()no mainargumento da plot()função.

Aqui está a parte relevante do meu código para plotagem:

# Define plotting parameters
pch_values <- c(3, 4, 18, 20)  # Symbols for different cases
colors <- c("red", "blue", "gray60")  # Colors for different noises
plot_titles <- c("Bias Curves", "Variance Curves")
coeff_labels <- c(expression(beta[0]), expression(beta[1]))


par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))  # 2x2 layout

for (plot_type in c("bias", "variance")) {
  for (coef_idx in seq_along(coeff_labels)) {
    coef_name <- ifelse(coef_idx == 1, "beta_0", "beta_1")
    
    # Compute y-axis range
    ylim_range <- range(sapply(results_coef, function(noise_list) {
      sapply(noise_list, function(case_list) {
        y_values <- abs(case_list[[plot_type]][[coef_name]])  # Absolute for bias
        return(y_values)
      })
    }), na.rm = TRUE)
    
    # Initialize the plot
    plot(NULL, log = "x", xlim = range(sample_sizes), ylim = ylim_range,
         xlab = "Sample Size (T)", ylab = str_to_title(plot_type),
         main = paste(plot_titles[ifelse(plot_type == "bias", 1, 2)], " (", coeff_labels[coef_idx], ") ", sep = ""))
    
    
    for (noise_idx in seq_along(noise_types)) {
      noise <- noise_types[noise_idx]
      case_counter <- 1
      
      for (scale in unique(unlist(lapply(cases, `[[`, noise)))) {
        y_values <- abs(results_coef[[noise]][[as.character(scale)]][[plot_type]][[coef_name]])  # Absolute for bias
        
        lines(sample_sizes, y_values, col = colors[noise_idx], lty = noise_idx, lwd = 1)
        points(sample_sizes, y_values, col = colors[noise_idx], pch = pch_values[case_counter])
        case_counter <- case_counter + 1
      }
    }
    
    # Add legends
    legend("topright", legend = str_to_title(noise_types), col = colors,
           lty = 1:length(noise_types), lwd = 1, bty = "n", inset = c(.2, 0))
    legend("topright", legend = paste("Case", seq_along(pch_values)),
           pch = pch_values, col = "black", bty = "n")
  }
}

Entretanto, em vez de exibir as expressões matemáticas \beta_0 e \beta_1 nos títulos dos gráficos, o maintexto mostra beta[0]e beta[1]como texto simples.

Pergunta:

Como posso formatar corretamente os títulos dos gráficos para incluir texto (por exemplo, "Curvas de polarização") e expressões matemáticas no estilo LaTeX como \beta_0​ e \beta_1​ no mainargumento de plot()?

Quaisquer sugestões ou abordagens alternativas seriam muito apreciadas.

Editar:

Aqui está um exemplo mínimo de trabalho:

# Define plotting parameters
coeff_labels <- c(expression(beta[0]), expression(beta[1]))
plot_titles <- c("Bias Curves", "Variance Curves")

# Attempting to plot with combined text and expressions
plot(NULL, xlim = c(1, 10), ylim = c(1, 10), 
     main = paste(plot_titles[1], " (", coeff_labels[1], ")", sep = ""))

Embora eu não tenha certeza se isso transmite corretamente o problema.

Em um experimento que estou realizando, a legenda do gráfico é criada da seguinte maneira:

eigenvalues <- t(c(1,2))
plot(1)
# Legend
legend_text <- cbind(
  paste("λ1 =", eigenvalues[, 1]),
  paste("λ2 =", eigenvalues[, 2])
)
legend("topright", legend = legend_text, text.col = "black", bty = "n", cex = 1)

Porém, em vez de "λ1 ="Eu gostaria de usar lambda[1], mas quando eu uso recebo o seguinte erro:

legend_text <- cbind(
  expression(lambda[1] == eigenvalues[, 1]),
  expression(lambda[2] == eigenvalues[, 2])
)

Error in cbind(expression(lambda[1] == eigenvalues[, 1]), expression(lambda[2] ==  : 
  cannot create a matrix from type 'expression'

Como usar expression(lambda[1]a legenda sem receber uma mensagem de erro?

Estou escrevendo um livro usando Quarto e, ao tentar gerar um enredo de três subfiguras, encontrei um problema ao adicionar uma sublegenda a cada subtrama respectivo. Este é o pedaço de código relevante:

```{r,echo=FALSE, fig.height=3, fig.width=9}
#| label: fig-plots
#| fig-cap: "Plots" 
#| fig-subcap:
#|   - "Plot 1"
#|   - "Plot 2" 
#|   - "Plot 3"
#| layout-ncol: 1

set.seed(5)
x_values <- seq(-4, 4, length.out = 1000)
a <- -sqrt(3)
b <- sqrt(3)

# Create the curve of the uniform distribution

par(mfrow = c(1,3), mar=c(5,1,5,1))
curve(dunif(x, min = a, max = b), from = -4, to = 4,
      n = 10000,
      col = "darkblue",
      lwd = 3,
      ylim = c(0, 0.45),
      ylab = '',
      xlab = 'E.K. = -1.2',
      main = expression("Uniform(min=-" * sqrt(3) * ",max=" * sqrt(3) * ")")
)
dunif_density <- function(x) dunif(x, min = a, max = b)
polygon(c(x_values, rev(x_values)), c(rep(0, length(x_values)), dunif_density(x_values)), col = 'darkblue', density = 20, border = NA)

# Create the normal distribution

curve(dnorm(x, mean = 0, sd = 1), from = -4, to = 4,
      n = 10000,
      col = 'darkblue',
      lwd = 3,
      ylim = c(0, 0.45),
      ylab = '',
      xlab = 'E.K. = 0',
      main = expression("Normal(" * mu * "=0, " * sigma * "=1)")
      )
dnorm_density <- function(x) dnorm(x, mean = 0, sd = 1)
polygon(c(x_values, rev(x_values)), c(rep(0, length(x_values)), dnorm_density(x_values)), col = 'darkblue', density = 20, border = NA)

# Create the Logistic Distribution

curve(dlogis(x, location = 0, scale = 0.55), from = -4, to = 4,
      n = 10000,
      col = 'darkblue',
      lwd = 3,
      ylim = c(0, 0.45),
      ylab = '',
      xlab = 'E.K. = 1.2',
      main = expression("Logistic(" * alpha * "=0, " * beta * "=0.55)")
      )
dlogis_density <- function(x) dlogis(x, location = 0, scale = 0.55)
polygon(c(x_values, rev(x_values)), c(rep(0, length(x_values)), dlogis_density(x_values)), col = 'darkblue', density = 20, border = NA)
```

E esta é a figura resultante:

Como você pode ver claramente, há apenas uma sublegenda, e não três, como esperado. Qual é a causa disso? E como posso consertar isso?