## ----echo = FALSE-------------------------------------------------------------
knitr::opts_chunk$set(
 fig.width = 8 ,
 fig.height = 12,
 fig.align ='center'
)

## -----------------------------------------------------------------------------
library("mvdpd")
library("ggplot2")
library("reshape2")
library(tidyr)
library(knitr)  

## -----------------------------------------------------------------------------
data("twocenturies50")
ordered_data <- twocenturies50[order(twocenturies50[,1]),]
n <- nrow(twocenturies50)
p <- ncol(twocenturies50)

## -----------------------------------------------------------------------------
best0 <- ordered_data[1:200,]
best0 <- best0/3600

## -----------------------------------------------------------------------------
best1 <- ordered_data[c(1:180, (n-19):n),]
best1 <- best1/3600

## -----------------------------------------------------------------------------
best2 <- ordered_data[1:200,]
set.seed(1234)
for (i in 1:200) {
  nn <- rbinom(n=1, size=10, prob=0.1)
  if (nn){
    pos <- sample(1:p, size=nn)
    best2[i,pos] <- runif(1, 0, 3)*best0[i,pos]
    best2[i,] <- sort(best2[i,]) 
  }
}
best2 <- best2/3600

## -----------------------------------------------------------------------------
min_val <- min(best0, best1, best2) 
max_val <- max(best0, best1, best2) 

heatplot <- function(mat) {
  df <- melt(mat)
  colnames(df) <- c("Event", "FinishingTime", "Hours")  
  df$FinishingTime <- as.numeric(as.factor(df$FinishingTime))
  
  ggplot(df, aes(FinishingTime, Event, fill = Hours)) +
    geom_tile() +
    scale_x_continuous(
      breaks = 1:10,
      labels = 1:10,
      expand = expansion(mult = c(0.01, 0.01))
    ) +
    scale_y_discrete(
      expand = expansion(mult = c(0.03, 0.03))
    ) +
    scale_fill_gradient2(
      low = "blue",
      mid = "white",
      high = "red",
      midpoint = (min_val + max_val) / 2,
      limits = c(min_val, max_val),
      breaks = round(seq(min_val, max_val, length.out = 5),0)
    ) +
    labs(
      x = "Top 10 Finishing Times",
      y = "Events",
      fill = "Hours"
    ) +
    theme_classic(base_size = 14)
}

## ----fig=TRUE-----------------------------------------------------------------
heatplot(best0)

## ----fig=TRUE-----------------------------------------------------------------
heatplot(best1)

## ----fig=TRUE-----------------------------------------------------------------
heatplot(best2)

## -----------------------------------------------------------------------------
perform.analysis <- function(x) { 
  # Maximum likelihood
  ml <- mvogammaML(x=x, method="multivariate",
    initial=list(delta=c(18, rep(0.5, 9)), lambda=5))
  # Maximum composite likelihood (minimum CDPD with beta=0)
  cml <- mvogammaML(x=x, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum composite DPD
  cdpd1 <- mvogammaDPD(x=x, beta=0.1, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd2 <- mvogammaDPD(x=x, beta=0.2, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd3 <- mvogammaDPD(x=x, beta=0.3, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd5 <- mvogammaDPD(x=x, beta=0.5, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum (multivariate) DPD
  mdpd1 <- mvogammaDPD(x=x, beta=0.1, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd2 <- mvogammaDPD(x=x, beta=0.2, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd3 <- mvogammaDPD(x=x, beta=0.3, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd5 <- mvogammaDPD(x=x, beta=0.5, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))

  results.multivariate <- cbind(
    c(ml$delta, ml$lambda),
    c(mdpd1$delta, mdpd1$lambda),
    c(mdpd2$delta, mdpd2$lambda),
    c(mdpd3$delta, mdpd3$lambda),
    c(mdpd5$delta, mdpd5$lambda)
  )

  results.composite <- cbind(
    c(cml$delta, cml$lambda),
    c(cdpd1$delta, cdpd1$lambda),
    c(cdpd2$delta, cdpd2$lambda),
    c(cdpd3$delta, cdpd3$lambda),
    c(cdpd5$delta, cdpd5$lambda)
  )

  colnames(results.multivariate) <- c("ML", "DPD(0.1)", "DPD(0.2)",
                                     "DPD(0.3)", "DPD(0.5)")
  colnames(results.composite) <- c("CML", "CDPD(0.1)", "CDPD(0.2)",
                                  "CDPD(0.3)", "CDPD(0.5)")
  rownames(results.multivariate) <- c(paste0("delta", 1:10), "lambda")
  rownames(results.composite) <- rownames(results.multivariate) 
  res <- list(MDPDE=results.multivariate, MCDPDE=results.composite)
  return(res)
}

## -----------------------------------------------------------------------------
results0 <- perform.analysis(best0)
results1 <- perform.analysis(best1)
results2 <- perform.analysis(best2)

## -----------------------------------------------------------------------------
casewiseCDPD <- 100*abs(results1$MCDPDE-results0$MCDPDE)/results0$MCDPDE
cellwiseCDPD <- 100*abs(results2$MCDPDE-results0$MCDPDE)/results0$MCDPDE

## -----------------------------------------------------------------------------
# Define the methods (X-axis categories) in the exact order from the image
methods <- c("ML", "CML", "DPD(0.1)", "DPD(0.2)", "DPD(0.3)", 
             "DPD(0.5)", "CDPD(0.1)", "CDPD(0.2)", "CDPD(0.3)", "CDPD(0.5)")

mybarplot <- function(X1, X2) {
  # Create a data frame with the extracted percentages
df <- data.frame(Method = factor(methods, levels = methods),
                 Row_1 = X1, Row_2 = X2)

# Convert data to long format for ggplot2
df_long <- pivot_longer(df, cols = c(Row_1, Row_2), 
                        names_to = "Parameter", values_to = "Percentage")

# Create the Plot
p <- ggplot(df_long, aes(x = Method, y = Percentage, fill = Parameter)) +
  
  # Add grouped bars with black outlines
  # Matching position_dodge width and bar width removes the inner gap
  geom_bar(
    stat = "identity",
    position = position_dodge(width = 0.8), # Changed from 0.85 
    color = "black",       
    width = 0.8,                            # Changed from 0.75 to match dodge width
    linewidth = 0.6
  ) +
  
  # Add data labels on top of the bars
  geom_text(
    aes(label = sprintf("%.0f", Percentage)),
    position = position_dodge(width = 0.8), # Must match the dodge width in geom_bar
    vjust = -0.5,
    size = 4,           
    color = "black"
  ) +
  
  # Map colors and use mathematical expressions for the legend labels (lambda and delta)
  scale_fill_manual(
    values = c("Row_1" = "#4C72B0", "Row_2" = "#DD8452"),
    labels = c(expression(lambda), expression(delta[1]))
  ) +
  
  # Configure Y-axis: force origin at 0 and set max limit slightly above data
  scale_y_continuous(
    limits = c(0, max(df_long$Percentage) * 1.15), 
    breaks = seq(0, 100, 20),
    expand = expansion(mult = c(0, 0)) 
  ) +
  
  # Axis and Legend Labels
  labs(
    x = "Estimation Method",
    y = "Absolute relative changes (in %)",
    fill = "Parameter"
  ) +
  
  # Theme Modifications
  theme_classic(base_size = 14) +
  
  theme(
    # Axis tick values bold and larger
    axis.text.x = element_text(angle = 45, hjust = 1, color = "black", size = 14, face = "bold"),
    axis.text.y = element_text(color = "black", size = 14, face = "bold"),
    
    # Axis titles
    axis.title = element_text(face = "bold", size = 16),
    
    axis.line = element_line(color = "black", linewidth = 0.7),
    axis.ticks = element_line(color = "black", linewidth = 0.7),
    
    # Legend formatting
    legend.position = c(0.88, 0.88),
    legend.background = element_rect(fill = "white", color = NA),
    legend.box.background = element_rect(color = "black", linewidth = 0.5),
    legend.title = element_text(face = "bold", size = 14),
    legend.text = element_text(size = 15), 
    
    # Plot margins
    plot.margin = margin(t = 15, r = 15, b = 10, l = 10)
  )

  # Display the plot
  print(p)
}

## ----fig=TRUE-----------------------------------------------------------------
mybarplot(casewiseCDPD[11,], casewiseCDPD[1,])

## ----fig=TRUE-----------------------------------------------------------------
mybarplot(cellwiseCDPD[11,], cellwiseCDPD[1,])

## -----------------------------------------------------------------------------
kable(cbind(results0$MDPDE, results0$MCDPDE))

## -----------------------------------------------------------------------------
kable(cbind(results1$MDPDE, results1$MCDPDE))

## -----------------------------------------------------------------------------
kable(cbind(results2$MDPDE, results2$MCDPDE))