# Stochasticke simulacne metody
# Radoslav Harman, KAMS FMFI UK
# Cvicenie 4 (Generovanie spojitych vektorov) - strucny zaznam prikazov pre R


rnorm.2d <- function(n, mu1=0, mu2=0, sd1=1, sd2=1, rho=0) {
# Generátor realizacii vektora (X1,X2) zo zdruzene normalneho rozdelenia
# mu1=E(X1), mu2=E(X2), sd1^2=D(X1), sd2^2=D(X2), rho=cor(X1,X2)

 X1 <- rnorm(n); X2 <- rnorm(n)
 Y1 <- sd1 * X1 + mu1
 Y2 <- sd2 * (rho * X1 + sqrt(1 - rho^2) * X2) + mu2

 return(rbind(Y1, Y2))
}

res <- rnorm.2d(5000); plot(res[1,], res[2,], pch=19, cex=0.1)
res <- rnorm.2d(5000, rho=0.5); plot(res[1,], res[2,], pch=19, cex=0.1)
res <- rnorm.2d(5000, rho=-0.9); plot(res[1,], res[2,], pch=19, cex=0.1)


rsimplex <- function(n, m) {
# Generátor n realizacii z rovnomerneho rozdelenia na m-rozmernom simplexe
# Poznamka: klasicky, nie pravdepodobnostny simplex

  res <- matrix(ncol=n, nrow=m)
  for(i in 1:n) {
    X <- runif(m); Y <- sort(X)
    Z <- c(Y[1], Y[2:m] - Y[1:(m-1)])
    res[,i] <- Z
  }

  return(res)
}

res <- rsimplex(10000, 2); plot(res[1,], res[2,], pch=19, cex=0.1, asp=1) # +m=3,10


rsphere <- function(n, m) {
# Generátor n realizacii z rovnomerného rozdelenia na povrchu m-rozmernej gule

  res <- matrix(ncol=n, nrow=m)
  for(i in 1:n) {
    X <- rnorm(m)
    res[,i] <- X / sqrt(sum(X^2))
  }

  return(res)
}

res <- rsphere(10000, 2); plot(res[1,], res[2,], pch=19, cex=0.1, asp=1) # +m=3,4,10


rball <- function(n, m) {
# Generátor n realizacii z rovnomerného rozdelenia vo vnutri m-rozmernej gule

  res <- matrix(ncol=n, nrow=m)
  for(i in 1:n) {
    X <- rnorm(m)
    res[,i] <- runif(1)^(1/m) * X / sqrt(sum(X^2))
  }

  return(res)
}

res <- rball(10000, 2); plot(res[1,], res[2,], pch=19, cex=0.1, asp=1) # +m=3,4,10


rtorus <- function(n, R, r) {
# Evkin generator n rovnomernych bodov na hranici torusu s polomermi R, r

  rphitor <- function(n, R, r) {
  # Generuje n hodnot s rozdelenim uhla phi pre torus s polomermi R, r
    res <- rep(0, n); k <- 0; OK <- FALSE
    while(!OK) {
      U <- 2 * pi * runif(1); V <- (R+r) * runif(1)
      if(abs(R + r * cos(U)) > V){
        k <- k + 1; if(k==n) OK <- TRUE
        res[k] <- U
      }
    }
    return(res)
  }

  phi <- rphitor(n, R, r)
  psi <- 2 * pi * runif(n)
  x <- cos(psi) * (R + r * cos(phi))
  y <- sin(psi) * (R + r * cos(phi))
  z <- r * sin(phi)

  return(rbind(x,y,z))
}

H <- matrix(rnorm(9), ncol=3); U <- eigen(H %*% t(H))$vectors
plot(t((U %*% rtorus(10000, 5, 1))[1:2,]), pch=19, cex=0.2, asp=1)