También es posible hacerlo de forma recursiva. El código en general
EDICIÓN:
Como te comenté hay algunas cosas que se puede hacer más eficiente, por ejemplo las funciones (p
, UI
, t
pueden mejorar, etc) deberian estar declaradas fuera de calcular()
para no volveraquí dejo el código con algunos cambios que comentaré después:
import itertools
import numpy as np
from scipy import optimize
from scipy import integrate
def p(x,b,n,y,s): return 1.0/(x*(b*n+y*np.log(x)-s*b*x))
def Ui(x,x1,b,n,y,s): return (integrate.fixed_quad(p, x, x1, args=(b,n,y,s), n=5)[0]-1.0)
def t(x,x1,b,n,y,s): return optimize.newton(Ui, x, args=(x1,b,n,y,s), maxiter=5000)
def ii(x,b,n,y,s): return ((y/b)*np.log(x)-s*x+n)
def di(x,b,n,y,s): return (s*x-(y/b))*(y*np.log(x)+n*b-s*b*x)
def p2(x,b,n,y,s): return (s*x-n)/(x*((y*np.log(x)+n*b-s*b*x)**2))
def pe(x,b,n,y,s): return (-np.log(x))/(x*((y*np.log(x)+n*b-s*b*x)**2))
def calcular(B):
LA= np.array(np.array(B)[:,0])
y= LA[0]
LB= np.array(B[1])[0]
b= LB[0]
x=0.99; x1=1.0
X=[1]
for _ in range(1,15): x1=t(x,x1, b, n, y, s); x=x1; X.append(x)
Xu = X[1:]
Y=X[0:2];C=X[1:3];D=X[2:4];F=X[3:5];G=X[4:6];H=X[5:7];J=X[6:8]
K=X[7:9];L=X[8:10];Q=X[9:11];W=X[10:12];E=X[11:13];R=X[12:14];T=X[13:15]
fi=[ii(x, b, n, y, s) for x in Xu]
iteorica=np.matrix([[fi[0]],[fi[2]],[fi[3]],[fi[4]],[fi[5]],[fi[6]],[fi[6]],[fi[7]],[fi[8]],[fi[9]],[fi[10]],[fi[11]],[fi[12]],[fi[13]]])
vdi=[di(x, b, n, y, s) for x in Xu]
def taobeta(N,M): return integrate.fixed_quad(p2, M, N, args=(b,n,y,s), n=10)[0]
valtaobeta= list(itertools.starmap(taobeta, [(Y),(C),(D),(F),(G),(H),(J),(K),(L),(Q),(W),(E),(R),(T)]))
sumtaobeta = np.cumsum(valtaobeta)
multib=[vdi[i]*sumtaobeta[i] for i in range(len(vdi))]
def taogamma(N,M): return integrate.fixed_quad(pe, M, N, args=(b,n,y,s), n=10)[0]
valtaogamma=list(itertools.starmap(taogamma, [(Y),(C),(D),(F),(G),(H),(J),(K),(L),(Q),(W),(E),(R),(T)]))
sumtaogamma=np.cumsum(valtaogamma)
g=[vdi[i]*sumtaogamma[i] for i in range(len(vdi))]
jacobiana = np.matrix(np.column_stack((multib,g)))
return(iteorica, jacobiana)
n=763.0
s=762.0
itecho=np.matrix([[3],[6],[25],[73],[222],[294],[258],[237],[191],[125],[69],[27],[11],[4]])
B = np.matrix([[0.5673],[0.0026]])
for i in range(3):
print('Iteración {}, valor de B:\n {}'.format(i, B).decode('UTF-8'))
iteorica, jacobiana = calcular(B)
B = (np.dot(np.linalg.inv(jacobiana.T*jacobiana),jacobiana.T*(iteorica-itecho)))+B
Los cambios más significativo son:
Limpiar un poco los import como te comenté, es especialmete importante evitar
from module import *
, sobretodo con módulos extensos y complejos comoNumPy
oSciPy
. Puedes terminar reescribiendo funciones de estos módulos sin querer y con resultados desagradables en general.En muchas de tus operaciones es plausible usar list-comprehensions evitando variables no necesarias y mejorando la eficiencia, por ejemplo:
fi=[] for x in Xu: qq=ii(x, b, n, y, s) fi.append(qq)
Puede quedar simplemente como:
fi=[ii(x, b, n, y, s) for x in Xu]
Haces varias sumas acumuladas usando ciclos
for
, puedes hacer esto eficientemente usandonumpy.cumsum
. Si usaras una versión de Python >= 3.2 podrias usar tambiénitertools.acumulate
. De esta forma puedes simplificar eficientemente cosas como:suma=0 sumtaobeta=[] for i in valtaobeta: suma=suma+i sumtaobeta.append(suma)
Por:
sumtaobeta = np.cumsum(valtaobeta)
Para calcular
jacobiana
puedes usarnp.column_stack
, reduciendo:jacobiana = matrix([[multib[0],g[0]],[multib[1],g[1]],[multib[2],g[2]],[multib[3],g[3]], [multib[4],g[4]],[multib[5],g[5]],[multib[6],g[6]],[multib[7],g[7]], [multib[8],g[8]],[multib[9],g[9]], [multib[10],g[10]],[multib[11],g[11]],[multib[12],g[12]],[multib[13],g[13]]])
a declararlas en cada iteraciónsimplemente:
jacobiana = np.matrix(np.column_stack((multib,g)))
- Todas las funciones que sean fáciles de sacar de
calcular()
están ahora fuera, solo quedan dentrotaobeta
ytaogamma
porque usas itertools.starmap y no es posible pasar los argumentosb
,n
,y
ys
a través de ella. Se podrian sacar también on algunos cambios pero no creo que merezca la pena. Esto evita que cada vez que llamemos a calcular de redefinan todas estas funciones, ahora son globales y no están dentro decalcular
. Nota que la segunda definición dep
que haces en tu código (la que usataobeta
) ahora se llamap2
.
Esto son solo ideas por si te interesa aplicarlas. No quiero modificar demasiadoSeguramente se pueden hacer otras mejoras pero hay conceptos del procedimiento matemático de tu código sin saber que hace cada cosa pero hay algunas cosaslógicamente se me escapan por lo que no termino de entender, por ejemplome atrevo a la función UI necesita 5 parámetros de entrada pero solo usas uno dentro de ella. Tampoco deberias importar NumPytocar más sin saber el significado de 4 formas distintas, esto termina siendo confuso y peligroso, sobretodo si usas la forma from module import *
cada variable/función.