Semana #3 y #4:

Jorge Rivera CC: 1031164050

Función para hallar el área bajo la curva de una integral y uso de librerias ggpubr y plotly.

library(agricolae)
library(ggplot2)
library(ggpubr)
library(plotly)
library(tidyr)
library(plyr)
library(readxl)

La integral a simular es la siguiente:

\[\int_{0}^{2e}ln(x)dx = \lim_{b \to 0^+}\int_{b}^{2e}ln(x)dx \]

El área real se halla mediante el cálculo del siguiente límite:

\[[xln(x)-x]_0^{2e}=2e*ln(2e)-2e\]

(area_exacta <- 2*exp(1)*log(2*exp(1)) - 2*exp(1))

## [1] 3.768339

Aproximaciones por medio de simulación:

Gráfico:

plot(log, xlim = c(0, 6), ylab = 'ln(x)')
abline(v = 0, lwd = 2) # asintota vertical
abline(v = 1, lty = 2, col = 'steelblue')
abline(h = 0, lwd = 2)
text(0.6, -1.5, 'Parte \n negativa', cex = 0.7)
text(4, 0.7, 'Parte \n positiva', cex = 0.7)

Simulación:

pn <- 100000
d1 <- NULL
dn1 <- NULL
x1 <- NULL
y1 <- NULL

for(i in 1:pn){
        x1[i] = runif(1, 0, 1)
        y1[i] = runif(1, -10, 0)
        d1[i] = y1[i] > log(x1[i])
        dn1[i] = ifelse(d1[i] == T, 1, 0)

}

plot(x1, y1, col = ifelse(dn1 == 1, 'steelblue', 'lightgray'), ylab = 'ln(x)', 
     main = 'Parte negativa', pch = 16, cex = 0.3)

prop1 <- length(dn1[dn1 == 1]) / pn # de todos los puntos del gráfico, muestreme los que son == 1 (azules)

pp <- 100000
d2 <- NULL
dn2 <- NULL
x2 <- NULL
y2 <- NULL

for(i in 1:pp){
        x2[i] = runif(1, 1, 2*exp(1))
        y2[i] = runif(1, 0, 2)
        d2[i] = y2[i] < log(x2[i])
        dn2[i] = ifelse(d2[i] == T, 1, 0)

}

plot(x2, y2, col = ifelse(dn2 == 1, 'steelblue', 'lightgray'), ylab = 'ln(x)', 
     main = 'Parte positiva', pch = 16, cex = 0.3, xlab = 'x')

prop2 <- length(dn2[dn2 == 1]) / pp # de todos los puntos del gráfico, muestreme los que son == 1 (azules)

Calculando área por simulación:

xn <- runif(pn, 1, 2*exp(1) )
yn <- runif(pn, 0, log(2*exp(1)))
xn1 <- runif(pn, 0, 1)
yn1 <- runif(pn, -2.9, 0)

a <- sum(xn > 1 & yn < log(xn))
b <- sum(xn1 < 1 & yn1 > log(xn1))

(area1 <- (log(2*exp(1)))*(2*exp(1) - 1)/pn * a)

## [1] 4.768162

(area2 <- (2.9)*(1)/(pn)*b)

## [1] 0.949866

(area_simulada <- area1 - area2)

## [1] 3.818296

Comparando:

area_exacta

## [1] 3.768339

area_simulada

## [1] 3.818296

Por medio de la simulación se obtuvo un valor aproximado de la integral, el cual puede mejorar su precisión con el aumento de puntos en la simulación.

---

---

Mejorando una funcion de la libreria agricolae agricolae:

mi_audps <- function(t, y){
  audps(y, t)
}

tiempos <- seq(7, 70, 7)
dano <- sort(runif(length(tiempos), 
                  min = 5, 80))
mi_audps(tiempos, dano)

## evaluation 
##   2407.711

mi_trat_audps <- function(yt){
  salida <- c()
  t <- unlist(yt[,1])
  for (i in seq(2,ncol(yt))) {
   salida[i] <-  mi_audps(t = t, y = unlist(yt[,i]))
  }
  salida <- salida[-1]
  return(salida)
}

df_da <- read_excel("~/Rprojects/UN/Computacion_e/datos_audps.xlsx")

audps_trat <- mi_trat_audps(df_da)

dfg <- data.frame(t = rep(df_da$tiempo, 3),
                  y = c(df_da$D_t1, df_da$D_t2, df_da$D_t3),
                  trat = rep(1:3, each=nrow(df_da)))
# dfg

Gráfico:

df_texto <- data.frame(label = round(audps_trat,2),
                       trat = 1:3)

ggplot(data = dfg, aes(x = t, y = y)) +
  geom_line() +
  geom_text(data = df_texto,
            mapping = aes(x = 20, y = 60, label = label)) +
  facet_wrap(~ trat)

Gráficos en plotly plotly:

data <- spread(dfg, trat, y)
data <- rename(data, c("1" = "Trat1", "2" = "Trat2", "3" = "Trat3"))

fig <- plot_ly(data, x = ~t, y = ~Trat1, type = 'scatter', mode = 'lines', name = 'Tratamiento 1')
fig <- fig %>% add_trace(y = ~Trat2, name = 'Tratamiento 2')
fig <- fig %>% add_trace(y = ~Trat3, name = 'Tratamiento 3')
fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9))

fig

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

p <- ggplot(dfg, aes(x = t, y = y, colour = as.factor(trat))) +
    geom_line() +
    facet_grid(. ~ trat)

fig <- ggplotly(p)
fig

Gráfico en ggpubr ggpubr:

dfg$trat <- factor(dfg$trat)
dfg$t <- rownames(dfg)
ggbarplot(dfg, x = 't', y = 'y',
          xlab = 'Tiempo', ylab = 'Daño',
          fill = "trat",               
          color = "black")