[Maxima-lang-es] [Fwd: Re: Ajuste no lineal por derivadas parciales de mínimos cuadrados]

[Argel G. A. <arg...@gm...>]

-------- Mensaje original --------
Asunto: 	Re: [Maxima-lang-es] Ajuste no lineal por derivadas parciales 
de mínimos cuadrados
Fecha: 	Mon, 13 Jul 2009 20:52:16 -0500
De: 	Argel Gastélum Arellánez <arg...@gm...>
A: 	Mario Rodriguez <bio...@te...>
Referencias: 	<4A5...@gm...> 
<4A5...@te...>



Mario Rodriguez escribió:
> Argel Gastélum Arellánez escribió:
>>     Hola, buenas tardes. Estoy tratando de hacer un ajuste no lineal 
>> de una función a un grupo de 214 parejas de datos, utilizando el 
>> método de derivadas parciales para encontrar mínimos cuadrados. Las 
>> instrucciones que estoy aplicando en wxMaxima 0.8.2 son las siguientes:
>>
>> --------------------------------------------------------------------------------- 
>>
>> M: read_matrix("C:/Users/Usuario/Desktop/Datos.txt")$
>> X: col(M,1)$
>> Y: col(M,2)$
>> Funcion_a_minimizar : '((1/length(X))*sum((F(X[i]) - Y[i])^2, i, 
>> 1,length(X)));
>> F(X) := exp(-a*X)/(b+c*X);
>> Suma_de_Cuadrados: first(''Funcion_a_minimizar);
>> Derivada_parcial_a: diff(Suma_de_Cuadrados,a)$
>> Derivada_parcial_b: diff(Suma_de_Cuadrados,b)$
>> Derivada_parcial_c: diff(Suma_de_Cuadrados,c)$
>> realonly: true$
>> Constantes_estimadas_Racionales: 
>> solve([Derivada_parcial_a=0,Derivada_parcial_b=0,Derivada_parcial_c=0],[a,b,c]); 
>>
>> Constantes_estimadas_Racionales2: 
>> (last(Constantes_estimadas_Racionales));
>> Constantes_estimadas_Decimales: 
>> float(last(Constantes_estimadas_Racionales));
>> [a[racional],b[racional],c[racional]]: map(rhs, 
>> Constantes_estimadas_Racionales2);
>> [a[decimal],b[decimal],c[decimal]]: map(rhs, 
>> Constantes_estimadas_Decimales);
>> load(draw)$
>> Escena_1:gr2d(
>> title="Curva ajustada a Datos Experimentales",
>> xlabel="X (unidades)",
>> ylabel="Y (unidades)",
>> xrange=[0,6],
>> yrange=[0,200],
>> key="Datos experimentales",
>> point_type=filled_circle,
>> point_size=1,
>> color=red,
>> points(M),
>> key="Curva Ajustada",
>> line_width=3,
>> color=blue,
>> explicit(ev(F(x),first(Constantes_estimadas_Racionales2)),x,0,10)
>> )$
>> wxdraw(Escena_1)$
>> /*
>> Ecuación de la Curva Ajustada:
>> */
>> y=exp(-a[racional]*x)/(b[racional]+c[racional]*x);
>> y=exp(-a[decimal]*x)/(b[decimal]+c[decimal]*x);
>> /*
>> Valores reportados por la fuente de los datos (Datos NIST Chwirut1):
>> a = 0.19027818370, b = 0.0061314004477, c = 0.010530908399.
>> */
>> --------------------------------------------------------------------------------- 
>>
>>
>>     Todo marcha bien hasta llegar a la orden:
>>
>> Constantes_estimadas_Racionales: 
>> solve([Derivada_parcial_a=0,Derivada_parcial_b=0,Derivada_parcial_c=0],[a,b,c]); 
>>
>>
>> donde se queda un rato trabajando el programa y después de unos 
>> segundos me arroja una serie de líneas como estas:
>> ----------------------------------------------------------------------------- 
>>
>> rat: replaced -1.0 by -1/1 = -1.0
>> rat: replaced -92.9 by -929/10 = -92.9
>> rat: replaced 0.5 by 1/2 = 0.5
>> ...
>> rat: replaced 6.0 by 6/1 = 6.0
>> Maxima encountered a Lisp error:
>>  Error in PROGN [or a callee]: The storage for BIGNUM is exhausted.
>> Currently, 7581 pages are allocated.
>> Use ALLOCATE to expand the space.
>> Automatically continuing.
>> To reenable the Lisp debugger set *debugger-hook* to nil.
>> ------------------------------------------------------------------------------ 
>>
>>     y a partir de aquí no avanza. He buscado pero no he encontrado la 
>> forma de expandir el espacio usando ALLOCATE y no tengo idea de qué 
>> se puede hacer. Esta secuencia de comandos me ha sido muy útil para 
>> ajustar curvas de Michaelis-Menten en experimentos de cinética 
>> enzimática y también para hacer regresiones lineales. He logrado 
>> encontrar los valores de "a", "b" y "c" usando el comando "lbfgs" (no 
>> me resulta muy práctico tener que dar valores iniciales, sólo por 
>> cuestiones de tiempo) y también con "lsquares" (que me ha parecido 
>> mucho mejor), pero me intriga el hecho de no poder hacerlo con este 
>> método de derivadas parciales.
>>
>>     Agradezco de antemano la ayuda brindada.
>>
>>     Saludos y buen inicio de semana.
>>
>> -- 
>>     Argel.
>>
>
> Hola Argel,
>
> La función 'solve' se utiliza para resolver ecuaciones y sistemas de 
> forma simbólica, lo que conlleva ciertas limitaciones prácticas y 
> teóricas.
>
> A simple vista no parece que tu sistema sea resoluble de forma 
> algebraica; incluso con tamaños muestrales mucho más pequeños te vas a 
> encontrar con polinomios de grado muy alto. Es el típico caso en el 
> que no hay más remedio que echar mano de procedimientos numéricos.
>
> Creo que tu mejor opción es utilizar las otras funciones a las que 
> haces referencia.
>
> Mario
>
>
   Hola Mario, muchas gracias por tu respuesta.

   ¿Entonces no es sólo un problema de falta de memoria disponible para 
wxMaxima?.

   Me provoca curiosidad porque comencé a usar estas instrucciones para 
encontrar mínimos cuadrados por derivadas parciales, debido a la 
necesidad de ajustar la función de Michaelis-Menten [V=Vmax*S/(Km+S)] a 
datos de velocidad de reacción contra concentración de sustrato (en este 
caso funcionó muy bien), y usando "lsquares_estimates" me salía un error 
como el siguiente:

-----------------------------------------------------------------------------------------------------
Coeficientes_de_la_curva:float(lsquares_estimates(M,[S,V],V=(Vmax*S)/(S+Km),[Vmax,Km]))
Division by 0
#0: lambda([x],ev(MSE,x))(x=[Vmax = 0,Km = -1/10])
#1: 
lsquares_estimates_exact(mse=('sum(('data[i,2]-'data[i,1]*Vmax/(Km+'data[i,1]))^2,i,1,10))/10,parameters=[Vmax,Km])
#2: 
lsquares_estimates(data=matrix([0.1,47.09],[0.5,130.3],[1,216.9],[1.5,263.1],[2,338.6],[3,311.5],[4,390.9],[5,428.8],[7,415....,variables=[S,V],equation=V 
= Vmax*S/(S+Km),parameters=[Vmax,Km],optional_args=[])(lsquares.mac line 
245)
-- an error.  To debug this try debugmode(true);
-----------------------------------------------------------------------------------------------------

   ¿En este caso qué se podría hacer para evitar el problema de la 
división por cero? ¿Será un problema del algoritmo utilizado en este 
método numérico?

   De nuevo muchas gracias por tu ayuda.

   Saludos y que estés bien.

--
   Argel.