First, let me explain two error messages that you got.
The first error
Error in -0.5 % * % t(mean_vec) % * % sigma_inv : non-conformable arguments
This error message shows that the cause of the error is the multiplication of arguments not conformable to the rule. The rule of the matrix multiplication A %*% B is that the number of column(s) of A if it's a matrix, or the length of A if it's a vector, has to equal the number of row(s) of B.
In density_fn the definition mean_vec<- matrix(c((x - mu_x), (y - mu_y))) makes mean_vec a matrix that has exactly 1 column and the number of rows equal to the total length of x combined with y. Consequently, if the number of rows of sigma_inv > 1,
the matrix multiplication t(mean_vec) %*% sigma_inv doesn't conform to the rule and will result in error. For example:
x <- 1:2 # length of 2
y <- 1:3 # length of 3
sigma <- matrix(1:4, nrow = 2)
sigma_inv <- solve(sigma)
mean_vec <- matrix(c(x - mean(x), y - mean(y)))
mean_vec # 1 row 5 columns
# [,1]
#[1,] -0.5
#[2,] 0.5
#[3,] -1.0
#[4,] 0.0
#[5,] 1.0
t(mean_vec) %*% sigma_inv
# Error in t(mean_vec) %*% sigma_inv : non-conformable arguments
It will not result in error if each of x and y has the length of 1. That's why this function does work when you use it for a single row. For example:
x <- 2 # length of 1
y <- 3 # length of 1
sigma <- matrix(1:4, nrow = 2)
mean_vec <- matrix(c(x - mean(x), y - mean(y)))
sigma_inv <- solve(sigma)
t(mean_vec) %*% sigma_inv
# [,1] [,2]
# [1,] 0 0
An alternative way to get the function to work properly for any length of x and y is to set the number of row and the number of column in mean_vec <- matrix() as follows:
mean_vec <- matrix(c((x - mu_x), (y - mu_y)),
ncol = nrow(sigma),
nrow = ncol(sigma))
and then to change the matrix multiplication -0.5 %*% mean_vec to a scalar multiplication -0.5 * mean_vec.
So, the function becomes:
density_fn <- function(x, y, mu_x, mu_y, sigma) {
mean_vec <- matrix(c((x - mu_x), (y - mu_y)),
ncol = nrow(sigma),
nrow = ncol(sigma))
sigma_det <- det(sigma)
sigma_inv <- solve(sigma)
frac <- 1 / (2 * pi * sqrt(sigma_det))
exponent <- exp(-0.5 * t(mean_vec) %*% sigma_inv %*% mean_vec)
den_fn <- frac * exponent
return(den_fn)
}
The second error
Error in match.fun(FUN) : c("'density_fn(flower$Width, flower$Length,
flower_w_mean, flower_l_mean, ' is not a function, character or
symbol", "' cov_matrix)' is not a function, character or symbol")
This error message shows that the value of FUN in apply is not properly specified. According to the documentation of apply, the value to be assigned to FUN should be
typically is either a function or a symbol (e.g., a backquoted name)
or a character string specifying a function to be searched for from
the environment of the call to apply
It means you should mention only the function name to FUN. The additional arguments of the function (the 2nd argument, the 3rd one,...) should be mentioned after FUN, not within FUN. Please check ?apply for details.
However, this step alone does not solve the problem because apply is suitable for a univariate input. Because your input is multivariate, mapply is more suitable. Other options include Map and loops using for etc.
How to apply functions using mapply
It is best explained by a simple example. Please check ?mapply for details. Suppose you have x and y and you want to get z = 2x + 3y. You want to vectorize the function to x and y.
x <- c(3,4,5)
y <- c(10,20, 30)
myfun <- function(x,y) 2*x + 3*y
z <- mapply(myfun, x, y)
z
#[1] 36 68 100
If you have arguments other than x and y in myfun, you should assign them to MoreArgs. Here is the example using density_fn.
mapply(density_fn, x = flower$Width, y = flower$Length, MoreArgs = list(
mu_x = flower_w_mean,
mu_y = flower_l_mean,
sigma = cov_matrix))
# [,1] [,2] [,3]
# [1,] 0.001293784 0.001000747 0.001929141
# [2,] 0.001293784 0.001000747 0.001929141
# [3,] 0.001293784 0.001000747 0.001929141
# [4,] 0.001293784 0.001000747 0.001929141
These steps do not produce errors. However, because of my limited knowledge on the subject, I do not guarantee that these steps properly represent the density function of bivariate normal distribution as you intended.
