R: búsqueda de raíces para un producto cartesiano de parámetros de función
Dada una función f (x, c, d) de x que también depende de algunos parámetros c y d. Me gustaría encontrar los ceros para un producto cartesiano de ciertos valores c_1, ..., c_n y d_1, ..., d_m de los parámetros, es decir, un x_ij tal que f (x_ij, c_i, d_j) = 0 para i = 1, ..., ny j = 1, ..., m. Aunque no es tan crucial, estoy aplicando un algoritmo de Newton-Raphson para la búsqueda de la raíz:
newton.raphson <- function(f, a, b, tol = 1e-5, n = 1000){
require(numDeriv) # Package for computing f'(x)
x0 <- a # Set start value to supplied lower bound
k <- n # Initialize for iteration results
# Check the upper and lower bounds to see if approximations result in 0
fa <- f(a)
if (fa == 0.0){
return(a)
}
fb <- f(b)
if (fb == 0.0) {
return(b)
}
for (i in 1:n) {
dx <- genD(func = f, x = x0)$D[1] # First-order derivative f'(x0)
x1 <- x0 - (f(x0) / dx) # Calculate next value x1
k[i] <- x1 # Store x1
# Once the difference between x0 and x1 becomes sufficiently small, output the results.
if (abs(x1 - x0) < tol) {
root.approx <- tail(k, n=1)
res <- list('root approximation' = root.approx, 'iterations' = k)
return(res)
}
# If Newton-Raphson has not yet reached convergence set x1 as x0 and continue
x0 <- x1
}
print('Too many iterations in method')
}
La función real que me interesa es más complicada, pero el siguiente ejemplo ilustra mi problema.
test.function <- function(x=1,c=1,d=1){
return(c*d-x)
}
Entonces para cualquier c_i y d_j puedo calcular fácilmente el cero por
newton.raphson(function(x) test.function(x,c=c_i,d=d_j),0,1)[1]
que aquí obviamente es solo el producto c_i * d_j. Ahora intenté definir una función que encuentre para dos vectores dados (c_1, ..., c_n) y (d_1, ..., d_m) los ceros para todas las combinaciones. Para esto, traté de definir
zeroes <- function(ci=1,dj=1){
x<-newton.raphson(function(x) test.function(x,c=ci,d=dj),0,1)[1]
return(as.numeric(x))
}
y luego use la función externa, por ejemplo,
outer(c(1,2),c(1,2,3),FUN=zeroes)
Desafortunadamente esto no funcionó. Recibí un mensaje de error
Error during wrapup: dims [product 6] do not match the length of object [1]
También podría haber una solución mucho mejor para mi problema. Estoy feliz por cualquier aportación.