# Distribution
Fx <- function(x) pnorm(x)
fx <- function(x) dnorm(x)


n <- 10
# ATTENTION: 6M simulation may take some time
z <- replicate(1e6, min(rnorm(n)))

# Probabiltiy z > y

prob <- function(x) {
  sapply(x, function(u) mean(z>u))
}


# Density
fos_k <- function(x, n, k) {
  fact <- factorial(n) / (factorial(k-1) * factorial(n-k))
  fact * fx(x) * Fx(x)^(k-1) *(1-Fx(x))^(n-k)
}

curve(fos_k(x, n=n, k=1), from=-1, to=1)
curve(prob, add=T, col=2, from=-1, to=1)