There is a surface defined by points in the form of a two-dimensional array. That is, the array indices are the x and y coordinates, and the values of the array elements are the z coordinate values of the corresponding points.
It is necessary to find a function a · x + b · y + c = z that is optimal from the point of least squares, that is, calculate the corresponding coefficients a, b, and c. Is there a way to do that in R?
Thank you!
This is the problem of finding the best fitting plane to a collection of points in 3D space.
There's a simple way to look at this problem. Consider each element in the matrix as an "observation" with two independent variables (x and y) and one dependent variable (z).
Then you can find your coefficients by least squares by simply running the linear model:
lm(z ~ x_indices + y_indices)
In fact, we can use this to construct a simple function that takes a matrix as input and gives the values of a, b and c as output:
best_plane <- function(any_2d_matrix)
{
x_values <- rep(seq(ncol(any_2d_matrix)), each = nrow(any_2d_matrix))
y_values <- rep(seq(nrow(any_2d_matrix)), ncol(any_2d_matrix))
z <- as.vector(any_2d_matrix)
suppressWarnings(result <- summary(lm(z ~ x_indices + y_indices))$coef)
return(c(a = result[2, 1], b = result[3, 1], c = result[1, 1]))
}
To see this at work, let's construct a matrix using pre-set values of a, b and c:
a <- 3.2
b <- 0.5
c <- -4
x_indices <- rep(1:10, each = 10)
y_indices <- rep(1:10, 10)
my_matrix <- matrix(a * x_indices + b * y_indices + c, nrow = 10)
my_matrix
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] -0.3 2.9 6.1 9.3 12.5 15.7 18.9 22.1 25.3 28.5
#> [2,] 0.2 3.4 6.6 9.8 13.0 16.2 19.4 22.6 25.8 29.0
#> [3,] 0.7 3.9 7.1 10.3 13.5 16.7 19.9 23.1 26.3 29.5
#> [4,] 1.2 4.4 7.6 10.8 14.0 17.2 20.4 23.6 26.8 30.0
#> [5,] 1.7 4.9 8.1 11.3 14.5 17.7 20.9 24.1 27.3 30.5
#> [6,] 2.2 5.4 8.6 11.8 15.0 18.2 21.4 24.6 27.8 31.0
#> [7,] 2.7 5.9 9.1 12.3 15.5 18.7 21.9 25.1 28.3 31.5
#> [8,] 3.2 6.4 9.6 12.8 16.0 19.2 22.4 25.6 28.8 32.0
#> [9,] 3.7 6.9 10.1 13.3 16.5 19.7 22.9 26.1 29.3 32.5
#> [10,] 4.2 7.4 10.6 13.8 17.0 20.2 23.4 26.6 29.8 33.0
Now let's see if we can retrieve our coefficients:
best_plane(my_matrix)
#> a b c
#> 3.2 0.5 -4.0
And if we add lots of random noise, we still come fairly close:
best_plane(my_matrix + rnorm(100))
#> a b c
#> 3.2486162 0.5054093 -4.3669805