"aniv" <- function(n, p){ if(missing(n) && missing(p)) error("um dos argumentos, n ou p deve ser fornecido") if(!missing(n) && !missing(p)) error("apenas um dos argumentos, n ou p deve ser fornecido") Prob <- function(n) 1 - exp(sum(log(365:(365-n+1))) - n*log(365)) VecProb <- Vectorize(Prob, "n") if(missing(n)) res <- sapply(p, function(y) which((VecProb(1:366) - y) > 0)[1]) if(missing(p)) res <- VecProb(n) return(res) } aniv(n=23) aniv(n=c(36,41, 86)) aniv(n=c(10, 20, 35, 50, 57)) aniv(n=366) plot(1:366, aniv(n=1:366), type="l", xlab="n", ylab="P[Coincidencia]") aniv(p=0.5) aniv(p=c(0.2, 0.4, 0.5, 0.7, 0.9, 0.99)) plot(1:100, aniv(n=1:100), type="l", xlab="n", ylab="P[Coincidencia]") arrows(c(1,aniv(p=0.5)),c(0.5, 0.5),c(aniv(p=0.5),aniv(p=0.5)),c(0.5,0), length=0.1) text(1, 0.5, 0.5, pos=2, off=0.1, cex=0.7) text(aniv(p=0.5),0 ,aniv(p=0.5), pos=1, off=0.2, cex=0.7) ## ##O problema das sequĂȘncias de caras e coroas ## "nTenta" <- function(N, padrao="HTT", media = TRUE){ padrao <- strsplit(padrao, NULL)[[1]] nc <- length(padrao) nTenta <- numeric(N) for(i in 1:N){ res <- sample(c("H","T"), nc, rep=T) n <- nc while(any(res != padrao)){ res <- c(res[2:nc], sample(c("H","T"), 1, rep=T)) n <- n+1 } nTenta[i] <- n } if(media) return(mean(nTenta)) else return(nTenta) } nTenta(100, "HTT", med=F) nTenta(27, "HTT", med=F) nTenta(100, "HTT") n1 <- nTenta(10000, "HTT", med=F) plot(table(n1)) nTenta(100000, "HTT") nTenta(100000, "HTH") ## ## O problema da carta premiada (Monty Hall) ## "jogo" <- function(){ cartas <- LETTERS[1:3] premio <- sample(cartas, 1) escolha <- sample(cartas, 1) sobra <- cartas[which(cartas != escolha)] mostra <- sample(sobra[which(sobra != premio)], 1) NTroca <- escolha Res.NT <- ifelse(NTroca == premio, "Ganhou", "Perdeu") Troca <- sobra[sobra != mostra] Res.T <- ifelse(Troca == premio, "Ganhou", "Perdeu") return(c(premio, escolha, mostra, NTroca, Res.NT, Troca, Res.T)) } jogo() set.seed(231) sim <- as.data.frame(t(replicate(10000, jogo()))) names(sim) <- c("premio", "escolha", "mostra", "NTroca", "Res.NT", "Troca", "Res.T") dim(sim) head(sim) prop.table(table(sim$Res.NT)) prop.table(table(sim$Res.T))