Estoy intentando realizar un ajuste global con el paquete symfit, siguiendo documentación de Symfit.

import numpy as np
import symfit as sf
import matplotlib.pyplot as plt
%matplotlib inline # for ipynb

# Generate example data
t = np.arange(0.0, 600.1, 30)
k = 0.005
C1_0, C2_0 = 1.0, 2.0
C1 = C1_0 * np.exp(-k*t)
C2 = C2_0 * np.exp(-k*t)

# Construct model
x_1, x_2, y_1, y_2 = sf.variables('x_1, x_2, y_1, y_2')
kg = sf.Parameter(value=0.01, min=0.0, max=0.1)
a_1, a_2 = sf.parameters('a_1, a_2')
globalmodel = sf.Model({
    y_1: a_1 * np.e**(- kg * x_1),
    y_2: a_2 * np.e**(- kg * x_2),
})

# Do fit
globalfit = sf.Fit(globalmodel, x_1=t, x_2=t, y_1=C1, y_2=C2)
globalfit_result = globalfit.execute()
print(globalfit_result)

### EDITED START
while globalfit_result.r_squared < 0.99:
    kg = sf.Parameter(value=globalfit_result.params['kg'])
    a_1 = sf.Parameter(value=globalfit_result.params['a_1'])
    a_2 = sf.Parameter(value=globalfit_result.params['a_2'])
    globalmodel = sf.Model({
        y_1: a_1 * np.e**(- kg * x_1),
        y_2: a_2 * np.e**(- kg * x_2),
    })
    globalfit = sf.Fit(globalmodel, x_1=t, x_2=t, y_1=C1, y_2=C2)
    globalfit_result = globalfit.execute()
### EDITED END

y_r = globalmodel(x_1=t, x_2=t, **globalfit_result.params)

# Plot fit
plt.plot(t,C1,'ro')
plt.plot(t,C2,'b+')
plt.plot(t,y_r[0],'r-')
plt.plot(t,y_r[1],'b-')
plt.show()

En este ejemplo, espero que el parámetro "kg" en el "modelo global" esté optimizado a 0,005. Sin embargo, el valor de "kg" es aproximadamente 9.6e-3, que es demasiado cercano al valor inicial (10.0e-3). Creo que hago algo estúpido, pero no puedo entenderlo.

¡Cualquier comentario y sugerencia es bienvenido!

EDITADO

Agregué un bucle while (muy feo) para obtener el mejor ajuste. No estoy seguro de por qué debería ser, pero parece funcionar.

1
heartfield 4 dic. 2016 a las 15:35

1 respuesta

La mejor respuesta

Parece que los límites están causando el problema. Los eliminé en mi prueba y luego todo funciona bien. Este es un problema conocido en symfit 0.3.3, y uno que ya están fijadas en la ̶ [̶ ̶m̶a̶s̶t̶e̶r̶ ̶] ̶ [1] ̶ rama sobre ̶G̶i̶t̶h̶u̶b̶.̶ ̶ ̶ i subido una nueva Dev versión que podría ahora instalar con ̶ ̶p̶i̶p̶ ̶i̶n̶s̶t̶a̶l̶l̶ ̶s̶y̶m̶f̶i̶t̶=̶=̶0̶.̶3̶.̶3̶.̶d̶e̶v̶1̶5̶5̶ ̶-̶-̶u̶p̶g̶r̶a̶d̶e̶ ̶, ̶ hasta que oficialmente liberación 0.3.4 ̶ (que será idéntica pero con Extended documentación) ̶, que ahora ha sido fijado en las versiones más nuevas.

Tenga en cuenta que cambié su np.e a sf.exp porque eso es simbólico. Mi código de trabajo está abajo, idéntico al suyo excepto por el cambio mencionado y ejecutándose en 0.3.3.dev155.

import numpy as np
import symfit as sf
import matplotlib.pyplot as plt

# Generate example data
t = np.arange(0.0, 600.1, 30)
k = 0.005
C1_0, C2_0 = 1.0, 2.0
C1 = C1_0 * np.exp(-k*t)
C2 = C2_0 * np.exp(-k*t)

# Construct model
x_1, x_2, y_1, y_2 = sf.variables('x_1, x_2, y_1, y_2')
kg = sf.Parameter(value=0.01, min=0.0, max=0.1)
a_1, a_2 = sf.parameters('a_1, a_2')
globalmodel = sf.Model({
    y_1: a_1 * sf.exp(- kg * x_1),
    y_2: a_2 * sf.exp(- kg * x_2),
})

# Do fit
globalfit = sf.Fit(globalmodel, x_1=t, x_2=t, y_1=C1, y_2=C2)
globalfit_result = globalfit.execute()
print(globalfit_result)

y_r = globalmodel(x_1=t, x_2=t, **globalfit_result.params)

# Plot fit
plt.plot(t,C1,'ro')
plt.plot(t,C2,'b+')
plt.plot(t,y_r[0],'r-')
plt.plot(t,y_r[1],'b-')
plt.show()
1
tBuLi 14 mar. 2017 a las 16:21
Gracias por su pronta respuesta con 0.3.3.dev155. Todo funciona bien excepto por el mensaje de error a continuación. C: \ Users \ User \ Miniconda2 \ lib \ site-packages \ scipy \ Optimize \ slsqp.py: 341: RuntimeWarning: valor no válido encontrado en mayor bnderr = donde (bnds [:, 0]> bnds [:, 1]) [0] C: \ Users \ User \ Miniconda2 \ lib \ site-packages \ symfit \ core \ fit.py: 1610: FutureWarning: la comparación con None dará como resultado una comparación de objetos por elementos en el futuro. si Ninguno en self.sigma_data.values ​​():
 – 
heartfield
5 dic. 2016 a las 05:37
Si lo entiendo correctamente, es una advertencia, ¿verdad? Entonces, ¿obtienes el resultado correcto? Me aseguraré de tener el segundo en 0.3.4;). La advertencia con respecto a los límites es más profunda y proviene del algoritmo de minimización, ¿todavía la recibe si, por ejemplo, cambia min a algo pequeño pero distinto de cero?
 – 
tBuLi
5 dic. 2016 a las 15:52
Si. ¡Obtuve el resultado correcto! Ambos son solo advertencias. La segunda advertencia permanece incluso si cambié el mínimo de "kg" a 0.001.
 – 
heartfield
6 dic. 2016 a las 05:57
Se deshizo de él y ahora se lanza la nueva versión que incluye documentos mejorados;)
 – 
tBuLi
6 dic. 2016 a las 19:12
1
PD. Si mi respuesta solucionó su problema, ¿le importaría aceptarla como la respuesta? :)
 – 
tBuLi
14 dic. 2016 a las 00:29