Em Rcpp
#include <RcppArmadillo.h>
#using namespace Rcpp;
#// [[Rcpp::depends(RcppArmadillo)]]
#// [[Rcpp::export]]
#List bitriang(int n = 100,double p = 0.5,double li = 0,
# double pm = 0.5,double ls = 1){
# double h1 = 2*p/(pm - li), h2 = 2*(1 - p)/(ls - pm);// Calculando as alturas
# int j = 0,i = 0;
# arma::vec ux, uy, I;
# // ux eh a uniforme do eixo x
# // uy eh a uniforme do eixo y
# // I eh a indicadora se foi aceito ou nao
# while(j<n){
# i = ux.n_elem;
# ux.resize(i+1);
# uy.resize(i+1);
# I.resize(i+1);
# ux(i) = (arma::randu<double>())*(ls - li) + li; // gera o escalar do eixo x
# uy(i) = (arma::randu<double>())*std::max(h1,h2);// gera o escalar do eixo y #que tambem eh o z
# //Primeira reta
# if(ux(i)<((li+pm)/2)){
# if(uy(i)<=((2*h1*(ux(i) - li))/(pm - li))){
# I(i) = 1;
# j = j + 1;
# }
# else{
# I(i) = 0;
# }
#
#}
# // segunda reta
# if((ux(i)>=((li+pm)/2))&&(ux(i)<pm)){
# if(uy(i)<=((2*h1*(ux(i) - pm))/(li - pm))){
# I(i) = 1;
# j = j + 1;
# }
# else{
# I(i) = 0;
# }
# }
# // terceira reta
# if((ux(i)>=pm)&&(ux(i)<((ls + pm)/2))){
# if(uy(i)<=((2*h2*(ux(i) - pm))/(ls - pm))){
# I(i) = 1;
# j = j + 1;
# }
# else{
# I(i) = 0;
# }
# }
# // quarta reta
# if((ux(i)>=((ls+pm)/2))&&(ux(i)<=ls)){
# if(uy(i)<=((2*h2*(ux(i) - ls))/(pm - ls))){
# I(i) = 1;
# j = j + 1;
# }
# else{
# I(i) = 0;
# }
# }
# i = i + 1;
# }
# return List::create(Named("z") = ux,
# Named("u") = uy,
# Named("I") = I,
# Named("h1") = h1,
# Named("h2") = h2);
#}
#// n = tamanho da amostra | por padrao n = 100
#// p = area do triangulo da esquerda | por padrao p = 0.5
#// li = limite inferior do eixo x | por padrao li = 0
#// pm = ponto medio do eixo x entre os triangulos | por padrao pm = 0.5
#// ls = limite superior do eixo x | por padrao ls = 1
#// --------------------------------------------------------------Teste: Cálculo das alturas dos triângulos.
#// H <- function(p,li,pm,ls){
#// return(c(h1 = (2*p)/(pm - li), h2 = 2*(1 - p)/(ls - pm)))
#// }
#// # Alturas para p = 0.7; li = 1; pm = 2; ls = 3
#// H(0.7,1,2,3)
#//
#//
#// plot.new()
#// plot.window(c(1,3), c(0,1.4))
#// axis(1, cex=.5)
#// axis(2, cex=.5)
#// mtext("f(x)", side=2, line=3, las=0, cex=1)
#// mtext("x", side=1, line=3, las=0, cex=1)
#//
#// curve(2*1.4*(x - 1),1,1.5,add = T)
#// curve(-2*1.4*(x - 2),1.5,2,add = T)
#// curve(2*0.6*(x - 2),2,2.5,add = T)
#// curve(-2*0.6*(x-3),2.5,3,add = T)
#// abline(h = 1.4,lty = 2);abline(h = 0.6,lty = 2)
#// abline(v = 1,lty = 2);abline(v = 2,lty = 2);abline(v = 3,lty = 2)Testando a Correlação e aderência da função.
#set.seed(4522)
#teste <- generateBitriAll(1000,0.7,1,2,3)
#hist(teste$z[teste$I==1],freq = FALSE,
#breaks = 30,xlab = "x",main = "")
#curve(2*1.4*(x - 1),1,1.5,add = T)
#curve(-2*1.4*(x - 2),1.5,2,add = T)
#curve(2*0.6*(x - 2),2,2.5,add = T)
#curve(-2*0.6*(x-3),2.5,3,add = T)
#abline(h = 1.4,lty = 2);abline(h = 0.6,lty = 2)
#abline(v = 1,lty = 2);abline(v = 2,lty = 2);abline(v = 3,lty = 2)A amostra é independente?
#par(mfrow = c(1,2))
#forecast::Acf(teste$z[teste$I==1],lag.max = 1000,
#xlab = "Defasagens",
#ylab = "Autocorrelação",
#main = "")
#forecast::Pacf(teste$z[teste$I==1],lag.max = 1000,
#xlab = "Defasagens",
#ylab = "Autocorrelação parcial",
#main = "")Gráfico que mostra as áreas que foram aceitas e rejeitadas, já que o método é de aceitação e rejeição.
#plot(teste$z,teste$u,
#col = ifelse(teste$I==1,"orange","lightblue"),
#pch = 20,xlab = "x",ylab = "f(x)")
#curve(2*1.4*(x - 1),1,1.5,add = T)
#curve(-2*1.4*(x - 2),1.5,2,add = T)
#curve(2*0.6*(x - 2),2,2.5,add = T)
#curve(-2*0.6*(x-3),2.5,3,add = T)
#abline(h = 1.4,lty = 2);abline(h = 0.6,lty = 2)
#abline(v = 1,lty = 2);abline(v = 2,lty = 2);abline(v = 3,lty = 2)