Plotting R function-CodePudding

I need to plot the number of students (Say 365) against the respective probability at least 2 are born on the same day.

I cannot seem to get it to plot the graph. Code runs but not sure.

nstudents <- function(x){

result = (1 - (factorial(365)/ 365^x*factorial(365-x)))
return(result)
}

x <- 2:365
y <- nstudents(x)

plot(x, y,
 main = "Number of Students VS Probability of atleast two born on same day",
 xlab = "x",
 ylab = "y",
 xlim = c(2, 365),
 ylim = c(0, 1),
 #type = "l",
 lwd = 5,
 col = "darkred")

This is the first bit of R code I have ever run so please do excuse my incompetence. Any assistance would be greatly appreciated.

Edit- After trying to print the output it appears it may be something to do with my formula. However I am not sure how to better formulate that as that is the exact formula.

Additional formula I am trying to use:

nstudents <- function(x){

result = (1 - exp(1)^((-x(x-1))/730))
print(result)
return(result)
}

However the above results in the following error:

"Error in x(x - 1): could not find function "x""

Is this because it is assuming x(x-1) is a function? Would I need to create a separate function for this?

CodePudding user response：

The problem lies in your function, which has a couple of issues. First of all, you have written:

result = (1 - (factorial(365)/ 365^x*factorial(365-x)))

Your demoninator should be 365^x * factorial(365-x), but the way you have written this means it is interpreted as (factorial(365) / 365^x) * factorial(365-x) because you have not explicitly wrapped the denominator in parentheses.

However, this doesn't really help, since factorial(365) is simply too large a number to represent in 64-bit floating point computers. You need to carry out the calculation in log format using lfactorial as MrFlick suggests, so your function would become:

nstudents <- function(n) {
  1 - exp(lfactorial(365) - n * log(365) - lfactorial(365 - n))
}

This allows:

x <- 2:365
y <- nstudents(x)

plot(x, y,
     main = paste("Number of Students VS Probability of at",
                  "least two born on same day"),
     xlab = "x",
     ylab = "y",
     xlim = c(0, 365),
     ylim = c(0, 1),
     col = "darkred")