Tengo 4 ecuaciones no lineales con tres incógnitas X, Y, y Z que quiero resolver. Las ecuaciones son de la forma:

F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ

…donde a, b y c son constantes que dependen de cada valor de F en las cuatro ecuaciones.

¿Cuál es la mejor manera de resolver esto?

  • Sólo para tu INFORMACIÓN: Es más común el uso de x, y y z para las variables independientes (es decir, los datos, en este caso), y a, b, c para los parámetros del modelo que usted está tratando de resolver. Cuando leí por primera vez su ecuación, estaba a punto de decir «pero eso lineal» (en términos de a, b y c). Sé que es una tontería discutir sobre la terminología, pero como es en la actualidad enunciado, muchas personas son propensas a errores de leer tu pregunta. (Bueno, la pregunta clara, sin embargo. +1)
  • También, es posible linealizar esta. Estoy escribiendo una respuesta, pero no tengo tiempo para terminarlo ahora. Si nadie contesta en el mientras tanto, voy a terminar mi respuesta, y publicarlo en una hora o dos (esperemos que alguien se me pegaba a él). Buena suerte!
  • El más perezoso manera (pero más fáciles de implementar, creo) es precompute para n (digamos 10) los valores para cada parámetro (para 1000 combinaciones en total), y de ver que la combinación de las puntuaciones más cercana a cero, y de zoom en alrededor de esa zona. Que deben trabajar con bastante facilidad por la mayoría de los tipos de ecuaciones, para darnos una idea de dónde buscar, pero hay más formas de fantasía que va a trabajar más rápido y(/o) sea más precisa.
  • scipy.optimize.brute hace exactamente lo que usted describe: docs.scipy.org/doc/scipy/reference/generated/…. Tenga en cuenta que usted necesita para buscar un 3D espacio de parámetros en este caso. Es muy sencillo, pero muy ineficiente. Dicho esto, si funciona, funciona. Si hay un montón de locales minimia y rangos de los parámetros son conocidos, puede ser un buen enfoque.
  • Cierto, pero en 3D es todavía muy fácil, y otro de los beneficios de la fuerza bruta que se obtiene una idea de la errorbars en su solución. (Lo que se dice, tan pronto como usted va sobre 3D, la fuerza bruta se convierte en desesperada)
  • Se puede utilizar sympy solver para hacer esto?

InformationsquelleAutor user1171835 | 2013-10-23

2 Comentarios

  1. 50

    Hay dos maneras de hacer esto.

    1. Uso no lineal solver
    2. Linealizar el problema y resolverlo en el de mínimos cuadrados sentido

    Instalación

    Así que, como yo entiendo su pregunta, usted sabe que F, a, b, y c en 4 puntos diferentes, y desea invertir para el modelo de parámetros X, y y Z. se tiene 3 incógnitas y 4 los datos observados puntos, así que el problema está sobredeterminada. Por lo tanto, vamos a estar resolviendo en el de mínimos cuadrados sentido.

    Es más común el uso de la frente de la terminología en este caso, así que vamos a voltear la ecuación alrededor. En lugar de:

    F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

    Vamos a escribir:

    F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

    Donde sabemos F, X, Y, y Z en 4 puntos diferentes (por ejemplo,F_0, F_1, ... F_i).

    Sólo estamos cambiando los nombres de las variables, no la ecuación de sí mismo. (Esto es más para mi la facilidad de pensar que cualquier otra cosa.)

    Lineal Solución

    Es posible linealizar esta ecuación. Usted puede resolver fácilmente para a^2, b^2, a b cos(c), y a b sin(c). Para hacer esto un poco más fácil, vamos a re-etiquetar las cosas de nuevo:

    d = a^2
    e = b^2
    f = a b cos(c)
    g = a b sin(c)

    Ahora la ecuación es mucho más simple: F_i = d + e X_i + f Y_i + g Z_i. Es fácil hacer un lineal de mínimos cuadrados de la inversión para d, e, f, y g. Podemos obtener a, b, y c de:

    a = sqrt(d)
    b = sqrt(e)
    c = arctan(g/f)

    Bueno, vamos a escribir esto en forma de matriz. Vamos a traducir 4 observaciones de (el código que voy a escribir va a tomar cualquier número de observaciones, pero vamos a mantenerlo de hormigón en el momento):

    F_i = d + e X_i + f Y_i + g Z_i

    En:

    |F_0|   |1, X_0, Y_0, Z_0|   |d|
    |F_1| = |1, X_1, Y_1, Z_1| * |e|
    |F_2|   |1, X_2, Y_2, Z_2|   |f|
    |F_3|   |1, X_3, Y_3, Z_3|   |g|

    O: F = G * m (soy un geophysist, así que nos G para «Funciones de Green» y m para «Parámetros del Modelo». Normalmente tendríamos que usar d para «datos» en lugar de F, así.)

    En python, esto se traduciría en:

    def invert(f, x, y, z):
        G = np.vstack([np.ones_like(x), x, y, z]).T
        m, _, _, _ = np.linalg.lstsq(G, f)
    
        d, e, f, g = m
        a = np.sqrt(d)
        b = np.sqrt(e)
        c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
        return a, b, c

    No-lineal de la Solución

    También se podría resolver esto mediante scipy.optimize, como @Joe sugirió. El más accesible de la función en scipy.optimize es scipy.optimize.curve_fit que utiliza una de Levenberg-Marquardt método por defecto.

    De Levenberg-Marquardt es una «escalada» algoritmo (bueno, va en descenso, en este caso, pero el término se usa de todos modos). En un sentido, hacer una estimación inicial de los parámetros del modelo (todos, por defecto en scipy.optimize) y siga la pendiente de observed - predicted en el espacio de parámetros de descenso a la parte inferior.

    Advertencia: Recoger el derecho de no-lineal de la inversión método de estimación inicial, y la optimización de los parámetros del método es muy «oscuro arte». Sólo se aprende haciendo, y hay un montón de situaciones en las que las cosas no funcionan correctamente. De Levenberg-Marquardt es un buen método si su espacio de parámetros es bastante suave (este debería ser). Hay un montón de otros (incluyendo los algoritmos genéticos, las redes neuronales, etc además de los métodos más comunes como el recocido simulado) que son mejores que en otras situaciones. No voy a ahondar en esa parte aquí.

    Hay un común gotcha que la optimización de los kits de herramientas de intentar corregir para que scipy.optimize no trate de manejar. Si su modelo de parámetros tienen diferentes magnitudes (por ejemplo,a=1, b=1000, c=1e-8), entonces usted necesitará cambiar la escala de las cosas por lo que son similares en magnitud. De lo contrario, scipy.optimize‘s «hill climbing» algoritmos (como la PELÍCULA) no calcular con precisión la estimación de la inclinación local, y le dará imprecisos resultados. Por ahora, estoy asumiendo que a, b, y c relativamente similares magnitudes. Además, ser conscientes de que, en esencia, todos los no-lineal métodos requieren para hacer una estimación inicial, y son sensibles a la que adivinar. Voy a dejar a continuación (acaba de pasar en como la p0 kwarg a curve_fit) porque el valor predeterminado a, b, c = 1, 1, 1 es bastante precisa de adivinar para a, b, c = 3, 2, 1.

    Con las salvedades de la forma, curve_fit espera que se pasa a una función, un conjunto de puntos en los que las observaciones fueron hechas (como una sola ndim x npoints de la matriz), y los valores observados.

    Así, si escribimos la función como esta:

    def func(x, y, z, a, b, c):
        f = (a**2
             + x * b**2
             + y * a * b * np.cos(c)
             + z * a * b * np.sin(c))
        return f

    Tendremos que envuelve a aceptar ligeramente diferentes argumentos antes de pasar a curve_fit.

    En pocas palabras:

    def nonlinear_invert(f, x, y, z):
        def wrapped_func(observation_points, a, b, c):
            x, y, z = observation_points
            return func(x, y, z, a, b, c)
    
        xdata = np.vstack([x, y, z])
        model, cov = opt.curve_fit(wrapped_func, xdata, f)
        return model

    Stand-alone Ejemplo de los dos métodos:

    Para darle una implementación completa, aquí un ejemplo en el que

    1. genera una distribución aleatoria de puntos para evaluar la función,
    2. evalúa la función en los puntos (usando el conjunto de parámetros del modelo),
    3. añade ruido a los resultados,
    4. y, a continuación, invierte para los parámetros del modelo utilizando los lineales y no lineales, métodos descritos anteriormente.

    import numpy as np
    import scipy.optimize as opt
    def main():
    nobservations = 4
    a, b, c = 3.0, 2.0, 1.0
    f, x, y, z = generate_data(nobservations, a, b, c)
    print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
    print linear_invert(f, x, y, z)
    print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
    print nonlinear_invert(f, x, y, z)
    def generate_data(nobservations, a, b, c, noise_level=0.01):
    x, y, z = np.random.random((3, nobservations))
    noise = noise_level * np.random.normal(0, noise_level, nobservations)
    f = func(x, y, z, a, b, c) + noise
    return f, x, y, z
    def func(x, y, z, a, b, c):
    f = (a**2
    + x * b**2
    + y * a * b * np.cos(c)
    + z * a * b * np.sin(c))
    return f
    def linear_invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)
    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c
    def nonlinear_invert(f, x, y, z):
    # "curve_fit" expects the function to take a slightly different form...
    def wrapped_func(observation_points, a, b, c):
    x, y, z = observation_points
    return func(x, y, z, a, b, c)
    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model
    main()
    • Esto es genial! Yo estaba buscando en el uso de scipy.optimizar el antes y no podía conseguir mi cabeza alrededor de ella. Sería de gran interés si no te importa tener un ir en él. Gracias de nuevo
    • Precioso! Este tipo de respuestas me recuerdan a los Stepanov de la cita de aquí: «érase una vez los programadores amaba las matemáticas y la conocía bien. (…) Hoy en día, tenemos programadores – incluso senior, director y jefe de programadores – que están orgullosos de no saber de matemáticas de la preparatoria. Se está poniendo de moda presumir de ser práctica, con las matemáticas sean vistos como académico de pacotilla. Creemos que la separación de la programación de las matemáticas es suicida para la programación. Matemáticamente personas analfabetas no innovar».
    • Gracias, me siento halagado! Excelente oferta!
  2. 2

    Usted probablemente querrá ser el uso de scipy no lineal de solucionadores de problemas, son muy fáciles: http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

    • Si usted no va más concreta de que, tal vez debería decir esto como un comentario, en lugar de una respuesta?
    • joe no tiene suficiente rep publicar comentarios todavía
    • uno de los puntos – que necesita para editar la respuesta, de lo contrario me puedo deshacer mi voto

Dejar respuesta

Please enter your comment!
Please enter your name here