Estimación de los valores de los parámetros de una ecuación diferencial

Question

Estoy tratando de ajustar el siguiente sistema ODE que describe la absorción y eliminación de fármacos.

Si $A_e$ es la cantidad de fármaco en estado de absorción y $A$ es la cantidad de droga en el cuerpo esto puede ser descrito como:

$$\begin{aligned} {dA_e \over dt} &= -k_a A_e\\ {dA \over dt} &= k_a A_e - k A\end{aligned}$$

Estoy tratando de estimar los valores de los parámetros la tasa de absorción ( $k_a$ ) y la tasa de eliminación ( $k$ ).

A continuación se muestra el código Matlab escrito para estimar los valores de los parámetros utilizando lsqnolin . Pero estoy recibiendo avisos de RelTol se ha aumentado a 2.2E-14 . ¿Puede alguien decirme qué es lo que falla en la forma en que lo he modelado?

Quiero saber como estimar usando un Sistema ODE (Es decir sin escribir la función objetivo como la solución de $A(t)$ ya que el modelo real tiene muchos otros estados)

time=[0.5;1;1.5;2;3;4;5;6;8;12];%time in hours

Drugdata=[0.33;1.23;1.95;2.72;3.51;3.63;3.47;3.22;2.39;1.13];

%parameters estimated:k_a and k
initial=[0.5,0.5];%initial guesses of the estimates

ODEinit=[8.48,Drugdata(1)];%initial values to the ODE. 

lb=[0,0];
ub=[100,100];
[fittedVals,resnorm,~,~,~,~,~]=lsqnonlin(@(xEstimate)errorFun(xEstimate,Drugdata,ODEinit,...
    time),initial,lb,ub);

fittedValues=model2(fittedVals',time,ODEinit);

close all
figure
h2=scatter(time,Drugdata,20,'b','filled');
hold on
h1=plot(time,fittedValues,'color', 'r');
ylabel('Drug')
xlabel('time')

function err=errorFun(xEstimate,DrugMeas,ODEinit,thours)
    Drug_est=model2(xEstimate,thours,ODEinit);%estiamted values of drug using the ODE model
    err=DrugMeas-Drug_est;
end

function Drugout= model2(xEstimate,thours,ODEinit)

options = odeset('AbsTol',1e-25,'RelTol',1e-25);
[~,values] = ode45(@(t,y)Equations(t,y),thours,ODEinit,options);

function s=Equations(~,y)
    %parameters estimated:k_a, k

    k=xEstimate(2);%rate of elimination
    ka=xEstimate(1);%rate of drug absorption

    s=zeros(2,1);
    s(1)=-ka*y(1);%absorption compartment
    s(2)=ka*y(1)-k*y(2);%drug

end
Drugout=values(:,2);
end

Answer 1

No ha sido la respuesta esperada. Era un método para ayudar a investigar el origen del problema.

Mientras tanto, Lutz Lehmann dio la pista que, de hecho, era una cuestión de especificación de la precisión de trabajo para el software.

No obstante, es interesante comparar los resultados de ambos métodos: Optimización de los parámetros del sistema de EDOs u optimización de los parámetros de la función que es solución de las EDOs.

Resolveremos analíticamente el sistema de EDO :

$${dA_e \over dt}=-k_a A_e\quad\implies\quad A_e=c\:e^{-k_a t}$$

$${dA \over dt}=k_a A_e-k\:A=k_a c\:e^{-k_a t}-k\:A$$ Resolución de la EDO lineal de primer orden $\quad {dA \over dt}+k\:A=k_a c\:e^{-k_a t}\quad$ conduce a : $$A(t)=c\:\frac{k_a}{k-k_a}e^{-k_a t}+Ce^{-k\,t}\tag 1$$ $c$ y $C$ son dos constantes.

Con $A(0)=0\quad\implies\quad C=c\:\frac{k_a}{k-k_a}$ $$A(t)=C\left(e^{-k\,t}-e^{-k_at} \right)\tag 2$$

Por lo tanto, esta función debe ajustarse aproximadamente a sus datos. Pruebe un programa de regresión para calcular los valores aproximados de $k$ , $k_a$ y $C$ .

Si este método falla al igual que el método de ajuste de ecuaciones diferenciales, se podría pensar que el modelo de EDOs no es del todo conveniente con respecto a los datos. Si tiene éxito se obtienen los valores aproximados de $k$ y $k_a$ que resuelve el problema.

EJEMPLO :

El ajuste de mínimos cuadrados medios para los cuatro parámetros de la ec. $(1)$ que es equivalente a : $$A(t)=c_1e^{-k_a t}-c_2e^{-k\,t}$$ da : $$k_a\simeq 0.251\quad;\quad k\simeq 0.335\quad;\quad c_1\simeq 37.961\quad;\quad c_2\simeq 39.450\quad;$$ RMSAE $\simeq 0.0879$

Los resultados de $k_a$ y $k$ son ligeramente diferentes a los resultados de Lutz Lehmann : $\quad k_a\simeq 0.252\quad;\quad k\simeq 0.301$

NOTA : El ajuste medio por mínimos cuadrados para los tres parámetros de la Ec. $(2)$ es menos bueno. Por lo tanto, suponiendo que $A(0)=0$ no era una buena idea.

Answer 2

Aparece una advertencia porque en el formato habitual de coma flotante de 64 bits una tolerancia relativa de 1e-25 no es factible. Redondear números reales a coma flotante ya da un error relativo de hasta 2e-16 a lo que hay que añadir la acumulación sobre los cálculos del método de solución y un error relativo de 2.2e-14 aparece como justificada para un caso mínimo para un valor positivo de la misma. Tus datos tienen 3 dígitos, por lo que una precisión de trabajo de 6 dígitos para la integración parece suficiente. Para ser más que exacto, establezca RelTol a 1e-9 y AbsTol a 1e-7 .

Haciendo experimentos con en la traducción del código a python usando scipy.optimize.least_squares, obtengo que

• Utilizando las tolerancias mínimas (es decir, las mayores) atol = 1e-4; rtol=1e-5; proporciona estimaciones de los parámetros [0.25206566 0.30065899]
• Uso de tolerancias moderadas atol = 1e-5; rtol=1e-7; proporciona estimaciones de los parámetros [0.25207333 0.30065716]
• Uso de tolerancias extremas atol = 1e-12; rtol=1e-13; proporciona estimaciones de los parámetros [0.25207395 0.30065723]

Como puede verse, los 3 o incluso 4 primeros dígitos de los parámetros estimados no dependen de las tolerancias de integración.

El gráfico del modelo ajustado es

Fijar la hora inicial en $0$ y haciendo que los valores iniciales sean también variables del proceso de optimización se obtiene el vector de parámetros ajustados $$ [k, k_a, A_{e0}, A_0] = [ 0.281374,\, 0.28137101,\, 10.95249421,\, -1.34117042] $$ con la trama

Forzar $A_0$ para ser positivo da un ajuste mucho peor.