a <- 2 # parametro de escala de la triangular

f.lenta <- function (x)
    ifelse (x <= -a | x >= a,
            0,
            ifelse (x > 0,
                    1/a^2 * (a - x),
                    1/a^2 * (a + x)))

f.rapida <- function (x)
    1/a^2 * ((a - x) * (x >= 0) * (x < a) +
             (a + x) * (x <  0) * (x > -a))

muestra.triangular.1 <- function (n=50, f=f.rapida)
{
    muestra <- numeric(n)
    i <- 1
    while (i <= n)
    {
        x <- runif (1, -a, a); y <- runif (1, 0, 1/a)
        if (y <= f(x))
        {
            muestra[i] <- x
            i <- i + 1
        }
    }
    muestra
}

muestra.triangular.2 <- function(n=50, f=f.rapida,
                                 muestra=numeric(0), crece=2)
    repeat {
        m <- crece * n # generamos de mas por los rechazos
        x <- runif(m, -a, a); y <- runif(m, 0, 1/a)
        muestra  <- c(muestra, x [y <= f(x)])
        if (length(muestra) >= n) return(muestra[1:n])
        m <- crece * m
    }

muestra.triangular.3 <- function(n=50)
{
    Y <- rbeta(n, 1, 2)
    V <- rbinom(n, 1, 0.5)
    a * 2 * Y * (V - 0.5)
}

muestra.triangular.4 <- function(n=50, f=f.rapida,
                                 I=8, graf=FALSE, extra=2) # I = num. interv.
{
    corte <- seq(-a, a, length.out = I + 1)
    altura <- sapply (1:I,
                      function(i)
                          optimize(f, corte[i:(i+1)], maximum=TRUE,
                                   tol=1e-6) $ objective)
    if (graf)
    {
        plot(f, -a, a, type = "l", col = 3)
        lines(c(-a,corte), c(0,altura,0), type = "s", col = 2)
    }
    anchura <- diff(corte)
    area <- sum(anchura * altura)
    m <- extra * n * area  # muestreamos de mas por los rechazos
    i <- sample(I, m, TRUE, altura)
    x <- anchura[i] * runif(m) + corte[i]
    y <- runif(m, 0, altura[i])
    x [y < f(x)] [1:n]
}

muestra.triangular.5 <- function(n=50, f=f.rapida,
                                 I=8, graf=FALSE, extra=2) # I = num. interv.
{
    corte <- seq(-a, a, length.out = I + 1) # equianchura supuesta
    altura <- sapply (1:I,
                      function(i)
                          optimize(f, corte[i:(i+1)], maximum=TRUE,
                                   tol=1e-6) $ objective)
    if (graf)
    {
        plot(f, -a, a, type = "l", col = 3)
        lines(c(-a,corte), c(0,altura,0), type = "s", col = 2)
    }
    anchura <- diff(corte)      # equianchura no supuesta
    C <- sum (anchura * altura) # integral de g(.)
    G <- c(0, cumsum (anchura * altura))
    G1 <- function (u)          # G^-1, inversa de funcion de distrib. G(.)
    {
        i <- sum(u > G)         # intervalo en que se encuentra "u"
        corte[i] + (corte[i+1]-corte[i])/(G[i+1]-G[i]) * (u-G[i])
    }
    m <- extra * n * C          # muestreamos de mas por los rechazos
    u <- runif (m, 0, C)
    A <- sapply(u, G1)
    g <- function(A) altura[ceiling((A-(-a))/anchura[1])] # equianchura supuesta
    b <- runif (m, 0, g(A))
    A [b < f(A)] [1:n]
}

system.time (hist (muestra.triangular.1(1e6, f=f.lenta), breaks=100, col=2))
system.time (hist (muestra.triangular.1(1e6),            breaks=100, col=3))
system.time (hist (muestra.triangular.2(1e6, f=f.lenta), breaks=100, col=4))
system.time (hist (muestra.triangular.2(1e6),            breaks=100, col=5))
system.time (hist (muestra.triangular.3(1e6),            breaks=100, col=6))
system.time (hist (muestra.triangular.4(1e6, f=f.lenta), breaks=100, col=7))
system.time (hist (muestra.triangular.4(1e6),            breaks=100, col=2))
system.time (hist (muestra.triangular.4(1e6, I=1),       breaks=100, col=3))
system.time (hist (muestra.triangular.4(1e6, I=100),     breaks=100, col=4))
system.time (hist (muestra.triangular.5(1e6),            breaks=100, col=5))
system.time (hist (muestra.triangular.5(1e6, I=1),       breaks=100, col=6))
system.time (hist (muestra.triangular.5(1e6, I=100),     breaks=100, col=7))
## > system.time (hist (muestra.triangular.1(1e6, f=f.lenta), breaks=100, col=2))
##    user  system elapsed 
##  26.536   0.008  27.384 
## > system.time (hist (muestra.triangular.1(1e6),            breaks=100, col=3))
##    user  system elapsed 
##  17.429   0.000  18.354 
## > system.time (hist (muestra.triangular.2(1e6, f=f.lenta), breaks=100, col=4))
##    user  system elapsed 
##    0.54    0.00    1.52 
## > system.time (hist (muestra.triangular.2(1e6),            breaks=100, col=5))
##    user  system elapsed 
##   0.591   0.000   1.420 
## > system.time (hist (muestra.triangular.3(1e6),            breaks=100, col=6))
##    user  system elapsed 
##   0.368   0.000   1.037 
## > system.time (hist (muestra.triangular.4(1e6, f=f.lenta), breaks=100, col=7))
##    user  system elapsed 
##   0.770   0.000   1.405 
## > system.time (hist (muestra.triangular.4(1e6),            breaks=100, col=2))
##    user  system elapsed 
##   0.516   0.000   1.108 
## > system.time (hist (muestra.triangular.4(1e6, I=1),       breaks=100, col=3))
##    user  system elapsed 
##   0.708   0.012   1.304 
## > system.time (hist (muestra.triangular.4(1e6, I=100),     breaks=100, col=4))
##    user  system elapsed 
##   0.463   0.000   1.004 
## > system.time (hist (muestra.triangular.5(1e6),            breaks=100, col=5))
##    user  system elapsed 
##   9.161   0.004   9.978 
## > system.time (hist (muestra.triangular.5(1e6, I=1),       breaks=100, col=6))
##    user  system elapsed 
##  15.032   0.052  16.028 
## > system.time (hist (muestra.triangular.5(1e6, I=100),     breaks=100, col=7))
##    user  system elapsed 
##   7.793   0.024   8.989