Esta es la función con miu = 0.5 y los lambdas iguales a 1

funcion_a_trozos <- function(x1, x2, b) {
  l1 <- 1
  l2 <- 1
  l3 <- 1
  miu <- 0.5
  
  if (x1 > (miu * x2)) {
    
    # si x1 > miu * x2
    
    fdp <- (b * l2 * x2^(b - 1)) * exp(-l2 * (x2^b)) * (b * (l1 + l3) * x1^(b - 1)) * exp(-((l1 + l3) * (x1^b)))
    # si x1 < miu * x2
    
  } else if (x1 < (miu * x2)) {
    fdp <- b * (l2 + l3 * (miu^b)) * (x2^(b - 1)) * exp(-((l2 + l3 * (miu^b)) * (x2^b))) * (b * l1 * (x1^(b - 1))) * exp(-(l1 * (x1^b))) 
    
     # si x1 = x2 = x
  } else if (x1 == (miu * x2)) {
    fdp <- ((l3 * (miu^b)) / (l1 * (miu^b) + l2 + l3 * (miu^b))) * b * (l1 + (l2 / (miu^b)) + l3) * (x1^(b - 1)) * exp(-((l1 + (l2 / (miu^b)) + l3) * (x1^b))) 
  }
  
  return(fdp)
}

GRÁFICOS A,C,E: son los gráficos para miu = 0.5 y betas 0.5, 1 y 5.

Este gráfico1 corresponde a beta= 0.5 , ver línea 98 gráfico 2 para beta = 1 . ver línea 135 gráfico 3 para beta = 5 , ver línea 169.

todos con miu=1

# gráfico a

library(rgl)

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos(x1[i], x2[j], 0.5)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig1<-persp(x1, x2, z, zlim=c(0,2), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,
      border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(a)") 

# SEGUNDO:

# gráfico c: para beta=1

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos(x1[i], x2[j], 1)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig2<-persp(x1, x2, z, zlim=c(0,1.7), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(c)")

# Tercero
#gráfico e

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos(x1[i], x2[j], 5)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig3<-persp(x1, x2, z, zlim=c(0,4.5), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(e)")

GRÁFICOS B,D,F

Creamos gráficos b,d,f son para miu = 1 y betas 0.5, 1 y 5:

funcion_a_trozos1 <- function(x1, x2, b) {
  l1 <- 1
  l2 <- 1
  l3 <- 1
  miu <- 1
  
  if (x1 > (miu * x2)) {
    
    # si x1 > miu * x2
    
    fdp <- (b * l2 * x2^(b - 1)) * exp(-l2 * (x2^b)) * (b * (l1 + l3) * x1^(b - 1)) * exp(-((l1 + l3) * (x1^b)))
    # si x1 < miu * x2
    
  } else if (x1 < (miu * x2)) {
    fdp <- b * (l2 + l3 * (miu^b)) * (x2^(b - 1)) * exp(-((l2 + l3 * (miu^b)) * (x2^b))) * (b * l1 * (x1^(b - 1))) * exp(-(l1 * (x1^b))) 
    
     # si x1 = x2 = x
  } else if (x1 == (miu * x2)) {
    fdp <- ((l3 * (miu^b)) / (l1 * (miu^b) + l2 + l3 * (miu^b))) * b * (l1 + (l2 / (miu^b)) + l3) * (x1^(b - 1)) * exp(-((l1 + (l2 / (miu^b)) + l3) * (x1^b))) 
  }
  
  return(fdp)
}
# gráfico b)

library(rgl)

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos1(x1[i], x2[j], 0.5)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig4<-persp(x1, x2, z, zlim=c(0,2), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,
      border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(b)") 

# SEGUNDO:

# gráfico d: para beta=1

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos1(x1[i], x2[j], 1)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig5<-persp(x1, x2, z, zlim=c(0,1.7), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(d)")

# Tercero
#gráfico f

# Define los valores de x1 y x2
x1 <- seq(0.1, 2.5, by = 0.05)
x2 <- seq(0.1, 2.5, by = 0.05)

# Crear una matriz para almacenar los valores de la función
z <- matrix(0, nrow = length(x1), ncol = length(x2))

# Calcular los valores de la función para cada par de (x1, x2)
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
    z[i, j] <- funcion_a_trozos1(x1[i], x2[j], 5)
  }
}

nrz=length(x1)
ncz=length(x2)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c( "lightblue") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

# Crear el gráfico de superficie 3D
fig6<-persp(x1, x2, z, zlim=c(0,4.5), theta = 45, phi = 5, expand = 0.6, col =color[facetcol], scale = TRUE, axes = TRUE,border = TRUE, box = TRUE, shade = 0.05, ticktype = "detailed", main = "(f)")

Funciones de densidad marginal:

marginal_x1_density<-function(x1,b){
l1 <- 1
l2 <- 1
l3 <- 1

fdm<-b*(l1+l3)*(x1^(b-1))*exp(-(l1+l3)*(x1^b))

return(fdm)
}
marginal_X2_density<-function(x2,b,u){
  l1<-1
  l2<-1
  l3<-1
  fdm<-b*(l2+l3*(u^b))*x2^(b-1)*exp(-(l2+l3*u^b)*(x2^b))
  return(fdm)
}
b <- c(0.5, 1, 1.5, 5)

# Crear el gráfico
plot(NULL, xlim = c(0, 2.5), ylim = c(0, 2.5), xlab = "x1", ylab = "f.d.m", main = "Función de densidad marginal (x1) para diferentes valores de beta", type = "n")

# Colores para las curvas
colores <- rainbow(length(b))

# Iterar sobre los valores de b y trazar cada curva
for (i in seq_along(b)) {
  x <- seq(0, 2.5, by = 0.01)
  y <- marginal_x1_density(x, b = b[i])
  
  # Añadir la curva al gráfico
  lines(x, y, col = colores[i], lty = i, lwd = 2)  # Usamos lty = i para que los tipos de línea varíen automáticamente
}

# Añadir la leyenda
legend("topright", legend = paste("b =", b), col = colores, lty = 1:length(b), lwd = 2, cex = 0.8)

b <- c(0.5, 1, 5)
u <- c(0.5, 1)

# Crear el gráfico
plot(NULL, xlim = c(0, 2.5), ylim = c(0, 2.3), xlab = "x2", ylab = "f.d.m", main = "Función de densidad marginal (x2) para diferentes valores de beta y u", type = "n")

# Colores para las curvas
colores <- rainbow(length(b) * length(u))


# Iterar sobre los valores de b y u y trazar cada curva
for (b_val in b) {
  for (miu in u) {
    x <- seq(0.1, 2.5, by = 0.05)
    y <- marginal_X2_density(x, b = b_val, u = miu)
    
# crear el dataframe
    
    
    # Añadir la curva al gráfico
    lines(x, y, col = colores[which(b == b_val) + length(b) * (which(u == miu) - 1)], lty = 1, lwd = 2)
  }
}

# Añadir la leyenda
legend("topright", legend = paste("b =", rep(b, each = length(u)), ", u =", rep(u, times = length(b))), col = colores, lty = 1, lwd = 2, cex = 0.8)

b <- c(0.5, 1, 1.5, 5)

# Colores para las curvas
colores <- rainbow(length(b))

df <- data.frame()

# Iterar sobre los valores de b y trazar cada curva
for (i in seq_along(b)) {
  x <- seq(0, 2.5, by = 0.01)
  y <- marginal_x1_density(x, b = b[i])
  
  df <- rbind(df, data.frame(x = x, y = y, beta = b[i]))  # Agregar beta al dataframe
  
}

# Crear el gráfico con ggplot fuera del bucle
ggplot(df, aes(x = x, y = y, color = factor(beta), linetype = factor(beta))) +
  geom_line(size = 1) +
  scale_color_manual(values = rainbow(length(b))) +
  scale_linetype_manual(values = 1:length(b)) +
  labs(x = "x1", y = "f.d.m", title = "Función de densidad marginal (x1) para diferentes valores de beta") +
  theme_classic() +
  theme(legend.position = "topright") +
  guides(color = guide_legend(title = "b"), linetype = FALSE)+
  coord_cartesian(ylim = c(0, 2.5))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

library(ggplot2)

b <- c(0.5, 1, 5)
u <- c(0.5, 1)

# Crear el dataframe para almacenar los datos
df <- data.frame()

# Iterar sobre los valores de b y u y trazar cada curva
for (b_val in b) {
  for (miu in u) {
    x <- seq(0.1, 2.5, by = 0.01)
    y <- marginal_X2_density(x, b = b_val, u = miu)
    
    # Agregar los datos al dataframe
    df <- rbind(df, data.frame(x = x, y = y, b = b_val, u = miu))
  }
}

# Crear el gráfico con ggplot
ggplot(df, aes(x = x, y = y, color = factor(b), linetype = factor(u))) +
  geom_line(size = 1) +
  scale_color_manual(values = colores, name = "b") +
  scale_linetype_manual(values = c(1,2), name = "u") +
  labs(x = "x2", y = "f.d.m", title = "Función de densidad marginal (x2) para diferentes valores de beta y u") +
  theme_classic() +
  theme(legend.position = "topright") +
  guides(color = guide_legend(title = "b"), linetype = guide_legend(title = "u"))