exercicios/lista 2

#1.Gerar uma permutação dos números 1; 2; 3; : : : ; n,
##considerando todas as n! possiveis permutações igualmente
##prováveis (exemplo 4b).

n = 10
Pi = c(1:n)

for(i in 1:n){
  ale = sample(i, 1:n, replace = FALSE)

  aux = round((ale)*runif(1)+1)
  if(aux != 1){
    aux2 = round((ale - 1)*runif(1)+1)
    Pi[aux] = Pi[aux2]
  }
};Pi

#2. Gerar valores de X ~ Geométrica(p) (exemplo 4b).
p = 0.4; n = 1000
inv = floor((log(runif(n))/(log(1-p))))
barplot(table(inv) +1)

#3. Implementar o algoritmo para gerar valores de X ~ Poisson(lamb) (seção 4.2).

val = runif(1)
n3 = 100

for(i in 1:n3){
  lamb = i
  pe = exp(-lamb)
  fe = pe
  if(val < fe){
    pe = (lamb*pe)/(i+1)
    fe = fe +pe
  }

}

# bernoulli

berinv = function( nc, p){
  X = c(1:nc)
  p1 = (1-p)

  for(i in 1:nc){
    y = runif(1)
    if(y < p1)
      X[i] = y
    else
      X[i] = 1
  }
  X
}
berinv(5,0.3)
#####################################################################################################################
# 2 - Cap.4 :
#1. Implementar algoritmo BINOMIAL da p. 56
#2. exercício 10
#3. exercício 12
#4. exercício 13
#5. exercício 14

######## inversa da binomial ##########

m = 1000
n = 15
p = 0.3

numbin = function(m, n, p){
  y = 0
  for (j in 1:m) {
    U = runif(n)
    x = 0

    for (i in 1:n) {

      if (U[i] <= 1 - p) {
        x[i] = 0
      } else {
        x[i] = 1
        }
    }
    y[j] = sum(x)
  }
  return(y)
}
test = numbin(m, n, p)
table(test)


######## 10.a ##########
##binomial negativa
#X ~ NBin (1; p) = Geo(p)
#X ~ NBin (r; p) = somatorio r i=1 Geo(p)
ng = 1000
r = 3
p = 0.5
tut = c()
geo = c()

for(j in 1:ng){

    tut[j] = sum (floor(((log(runif(r)))/(log(1-p)))) + 1)  # rgeom(1, p) ((log(runif(1)))/(log(p))) + 1

  #geo[j] = sum(tut)
}
barplot((table(tut)))


######## 10.b teorico ######

#p = j
#j = r, r+1, ...

p1 = 0.3
r1 = 5
vl = 0
for (i in 1:100) {
  j1 = r1 + i
  vl[i] = ((j1*(1-p1))/(j1+1-r1))*p1

}
((j(1-p))/(j+1-r))*p

r2 = 5
vl2 = 0
p2 = 0.3
for(i in 0:r2){
  j2 = r2 + i
  #p2 = j2
  vl2[i] = (factorial(j2))/((factorial(j2+1+r2)*factorial(r2-1)))*(p2^r2)*(1-p2)^(j2+1-r2)

}

# pj = (factorial(j))/((factorial(j+1+r)*factorial(r-1)))*(p^r)*(1-p)^(j+1-r)

n = 100
r = 5
j = r
p = 0.3
vla = 0
ja = c()
#pj = ((factorial(j-1))/((factorial(j-r)*factorial(r-1))))*(p^r)*(1-p)^(j+1-r)

for (i in 0:n) {
  #ja <- c()
  j = j+1
  ja[i] = j

  vla[i] = ((ja[i]*(1-p)*(ja[i]+1-r))/(ja[i]+1-r))*((factorial(ja[i]-1))/((factorial(ja[i]-r)*factorial(r-1))))*(p^r)*(1-p)^(ja[i]+1-r)


  }
plot(vla)


######## 10.c 4.3 ######

####### 12 #######
#prob condicional x>k

#gerando as probs
lamb = 3
k = 11
#x<k
#x = i
probabilidades = function(lamb, k){
  ia = 0
  for( i in 1:k-1){
    ia[i] = (exp(-lamb)*lamb^i)/factorial(i)
  }

  ja = 0
  for( j in 1:(k)){
    ja[j] = (((exp(-lamb)*lamb^j)/factorial(j)))
  }

  ij =0
  for (l in 1:k-1){
    ij[l] = ia[l]/sum(ja)
  }
  return( ij)
};pipa = probabilidades(lamb, k)
plot(pipa)
#
####gerando os num.
####seria pela inversa
X = 0
i = 0
j = 1

while (j < k){
  ru = runif(1)
  print(ru)
  if(ru < pipa[j]){
    X[j] = i
    j = j +1
  }else{
    i = i +1}
}

table(X)


####### 13 #####

####### 14 ####
tam = 100000

j = 0
X = 0
prob.vetor = c(0.11, 0.20, 0.31, 0.40, 0.51, 0.60, 0.71, 0.80, 0.91)
valor.vetor = 5:14

for(i in 1:tam){
    u = runif(1)
    if( u <= prob.vetor[1]){
      X[i] = valor.vetor[1]
    }else{
      if( u <= prob.vetor[2]){
        X[i] = valor.vetor[2]
      }else{
        if( u <= prob.vetor[3]){
          X[i] = valor.vetor[3]
        }else{
          if( u <= prob.vetor[4]){
            X[i] = valor.vetor[4]
          }else{
            if( u <= prob.vetor[5]){
              X[i] = valor.vetor[5]
            }else{
              if( u <= prob.vetor[6]){
                X[i] = valor.vetor[6]
              }else{
                if( u <= prob.vetor[7]){
                  X[i] = valor.vetor[7]
                }else{
                  if( u <= prob.vetor[8]){
                    X[i] = valor.vetor[8]
                  }else{
                    if( u <= prob.vetor[9]){
                      X[i] = valor.vetor[9]
                    }else{
                      X[i] = valor.vetor[10]
                    }
                  }
                }
              }
            }
          }
        }
      }
    }

  }

barplot(table(X)/tam)



# Lista 3 - Cap.5:
#1. Implementar algoritmo do exemplo 5f.
#2. exercício 1
#3. exercício 2
#4. exercício 6
#5. exercício 8a

#### gerar rv de uma normal ####

t = 1000
i = 1
lambida = 1


X = 0

while(i <= t){
  Y = rexp(1, lambida)
  U = runif(1)

  if(U <= exp((-(Y - 1)^2)/2)){
    X[i] = Y
    i = i +1
  }
}
h = round(X,2)
hist(h, freq = FALSE)

#### 01 ####

fx = function(x) exp(x)/(exp(1) - 1)
plot(fx)

u = runif(1)

t = log(u*(exp(1) - 1)+1)

m = 100000
oct = 0

for( i in 1:m){
  u = runif(1)
  oct[i] = log(u*(exp(1) - 1)+1)
}
hist(table(round(oct, 2))/m, freq = FALSE)
mean(oct)
var(oct)