rtriangle = function(n, a=0, b=1, c=0.5) {
# Use the inverse CDF method to generate observations from a triangular dist
# The parameters for the trianglulr distribution are the two endpoints
# of the triangle (a and b) and the location of the peak (c)
# n is the number of observations
# The output is a vector of n pseudo-random values from the triangular 
# distribution.

  u = runif(n)
# Determine the length of the base of the triangle and 
# the height of the density
  base = b-a
  Fc = (c-a)/base

  leftside = (u < Fc)
  t = u
# Find the IF(u) for those values of u with IF(u) between a and c and then
# find IF(u) for those values of u with IF(u) between c and b
  t[leftside] = a + sqrt(u[leftside]*base*(c-a))
  t[!leftside] = b - sqrt((1-u[!leftside])*base*(b-c))
  invisible(t)
}

dtriangle = function(x, a=0, b=1, c=0.5) {
# The parameters for the triangular distribution are the two endpoints
# of the triangle (a and b) and the location of the peak (c)
# x is the values for which to compute f(x), where f is the density function
# for the triangular distribution

  dx = rep(0, length(x) )
  dx[a < x & x <= c] = 2*(x[a<x & x<=c] - a)/
                                       ((b-a)*(c-a))
  dx[ c<x & x<b ] = 2*(b - x[c<x & x<b])/((b-a)*(b-c))
  dx
}


ARTri = function(n, a = 0, b = 1, c = 0.5, scale = 1, s1=4, s2 = 3)
{

# This function uses Acceptance rejection method to sample from a beta
# distribution.  
# Pairs of points (x,y) are generated uniformly over a triangle.
# The return value are those x of the (x,y) pairs where y is less than f(x)
# Here f is the density for the beta. 

# The inputs are the endpoints of the triangle (a,b) and the location of 
# the vertex c. (Note the triangle has base along the x-axis).
# The scale argument blows up the triangle to make it taller.
# Also the two parameters of the beta are given s1 and s2
# The value n is the number of observations to generate uniformly over the
# triangle.  Note that fewer than n values are returned when (x,y) pairs
# are rejected.


  x = rtriangle(n, a, b, c)
  y = runif(n, min=0, max = scale*dtriangle(x, a, b, c))
  keep = y < dbeta(x , s1, s2)
  plotCurves(a, b, c, scale, s1, s2)
  points(x[keep],y[keep],cex=0.5,col="blue")
  points(x[!keep],y[!keep],cex=0.5,col="red")
  invisible(x[keep])
}
  
plotCurves = function(a, b, c, scale = 1, s1=4, s2=3) {
# This function plots the triangle and the beta distribution.
# It is helpful to determine whether or not the triangle contains
# the beta distribution.
# See the function ARTri for a description of the input parameters.

  begin = min(a, 0)
  end = max(b, 1) 
  domain = seq(begin, end, (end-begin)/1000)
  betaCurve = dbeta(domain, s1, s2)
  triCurve = dtriangle(domain, a, b, c)
  ht = max(betaCurve, scale*triCurve)
  plot(domain, betaCurve, ylim=c(0, ht),
          type ="l", lwd= 2)
  points(domain, triCurve * scale, 
          type ="l", lwd =2)
}