#########################
A <- function(p)
{
   # Alias tabulka pre pravdepodobnostny vektor p
   # Hodnoty, ktore sa nadobudaju, su 1,2,...,length(p)
   
   n <- length(p); pup <- p
   blocked <- rep(FALSE, n)
   A <- matrix(0, nrow = n - 1, ncol = 4)
   
   for (k in 1:(n - 1)) {
      is <- which.min(pup + 2*blocked); blocked[is] <- TRUE
      js <- which.max(pup - 2*blocked)
      pdn <- rep(0, n); pdn[is] <- pup[is]*(n - k); pdn[js] <- 1 - pdn[is]
      A[k, 1] <- is; A[k, 2] <- pdn[is]; A[k, 3] <- js; A[k, 4] <- pdn[js]
      pup <- pup*(n - k)/(n - k - 1) - pdn/(n - k - 1)
   }
   
   A
}

r.alias <- function(N, p) {
   A <- A(p); V <- A[, c(1, 3)]; P <- A[, c(2, 4)]; n <- length(p)
   res <- rep(0, N); I <- ceiling((n - 1)*runif(N)) 
   for (i in 1:N) res[i] <- V[I[i], ceiling(runif(1) + 1 - P[I[i]])]
   res
}

######################
function(N,mu,x) 
{
# Odhad pravdepodobnosti udalosti, ze proces Sn niekedy prekroci x
# Sn=X1+...+Xn, Xi su iid N(-mu,1)

z<-rep(0,N)
for(i in 1:N){
 OK<-FALSE; S<-0
  while(!OK){
  S<-S+rnorm(1,mean=mu,sd=1)
  if(S>x){OK<-TRUE;z[i]<-exp(-2*mu*S)}
  }
 }
est<-mean(z)
del<-qnorm(0.975)*sd(z)/sqrt(N)
list(est-del,est+del)
}


####################
function(n,x,th,lam) 
{
# Odhad pravdepodobnosti exremalnej hodnoty compound Poisson sum S=X1+...+XN
# n ... pocet simulacii
# x ... prahova hodnota urcujuca rare event
# th ... parameter exponencialnej zmeny miery (od 0 do 1)
# lam ... parameter Poissonovho rozdelenia premennej N
# (Kazde Xi ma exponencialne rozdelenie s parametrom 1)

N<-rpois(n,lambda=lam)
res<-rep(0,n)
for(i in 1:n){
  S<-sum(rexp(N[i],rate=1-th))
  if(S>x) res[i]<-exp(-th*S)/(1-th)^N[i]
  }
print(summary(res))
est<-mean(res); delI<-qnorm(0.975)*sd(res)/sqrt(n)
list(P=est,Ilow=est-delI,Iup=est+delI)
}