Weibull distribution

x = seq(3, 3.5, by = 0.1)
matplot(x, cbind(((x-1)/x)^(1/x), factorial(1/x), (log(2))^(1/x)), type = "l", ylab = "")
legend("topleft",, c("Mode", "Mean", "Median"), col=1:3, lty=1:3)
f = function(k) ((k-1)/k)^(1/k) - factorial(1/k)
g = function(k) (log(2))^(1/k) - factorial(1/k)
h = function(k) ((k-1)/k)^(1/k) - (log(2))^(1/k)
ps = c(uniroot(f, c(3, 3.5))$root, uniroot(g, c(3, 3.5))$root,
  uniroot(h, c(3, 3.5))$root)
abline(v = ps, lwd = 2, lty = 2, col = "grey")
s = function(k) factorial(3/k)-3*factorial(1/k)*(factorial(2/k)-(factorial(1/k))^2) - (factorial(1/k))^3
uniroot(s, c(2, 50))