Si he entendido bién el problema deseas dando un valor inicial de B
recalcularla x veces. Para calcular B
usas la propia B
, jacobiana
, itecho
e iteorica
. Lo que puedes hacer es encerrar el código en una función que tenga como parámetro de entrada B
y que retorne iteorica y jacobiana, calculas B
y luego la vuelves a llamar con el nuevo valor de B
, algo así:
import scipy.optimize
from scipy import optimize
import numpy as np
from itertools import starmap
import itertools
from numpy import matrix
from numpy import *
from numpy.linalg import *
n=763.0; s=762.0
itecho=matrix([[3],[6],[25],[73],[222],[294],[258],[237],[191],[125],[69],[27],[11],[4]])
def calcular(B):
LA= np.array(np.array(B)[:,0])
y= np.array(LA[0])
LB= np.array(B[1])[0]
b= np.array(LB[0])
def p(x): return 1.0/(x*(b*n+y*np.log(x)-s*b*x))
def Ui(x,b,n,y,s): return (scipy.integrate.fixed_quad(p, x, x1, args=(), n=5)[0]-1.0)
def t(x, b, n, y, s): return scipy.optimize.newton(Ui, x, args=(b,n,y,s), maxiter=5000)
x=0.99; x1=1.0
X=[1]
Xu=[]
i=1
for v in range(1,15):
x1=t(x, b, n, y, s) ; x=x1; i=i+1; X.append(x); Xu.append(x)
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]
def ii(x,b,n,y,s): return ((y/b)*np.log(x)-s*x+n)
fi=[]
for x in Xu:
qq=ii(x, b, n, y, s); fi.append(qq)
iteorica=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]]])
def di(x,b,n,y,s): return (s*x-(y/b))*(y*np.log(x)+n*b-s*b*x)
vdi=[]
for x in Xu:
rr=di(x, b, n, y, s); vdi.append(rr)
def p(x,b,n,y,s): return (s*x-n)/(x*((y*np.log(x)+n*b-s*b*x)**2))
def taobeta(N,M): return scipy.integrate.fixed_quad(p, 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)]))
suma=0
sumtaobeta=[]
for i in valtaobeta:
suma=suma+i; sumtaobeta.append(suma)
multib=[]
for i in range(len(vdi)):
mul=vdi[i]*sumtaobeta[i]; multib.append(mul)
def pe(x,b,n,y,s): return (-np.log(x))/(x*((y*np.log(x)+n*b-s*b*x)**2))
def taogamma(N,M): return scipy.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)]))
suma=0
sumtaogamma=[]
for i in valtaogamma:
suma=suma+i; sumtaogamma.append(suma)
g=[]
for i in range(len(vdi)):
mul=vdi[i]*sumtaogamma[i]; g.append(mul)
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]]])
return(iteorica, jacobiana)
B = np.matrix([[0.5673],[0.0026]])
for i in range(3):
print('Iteración {}, valor de B:\n {}'.format(i, B))
iteorica, jacobiana = calcular(B)
B = (dot(inv(jacobiana.T*jacobiana),jacobiana.T*(iteorica-itecho)))+B
También es posible hacerlo de forma recursiva.
EDICIÓN:
Como te comenté hay algunas cosas que se pueden mejorar, aquí 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 como NumPy
o SciPy
. 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 usando numpy.cumsum
. Si usaras una versión de Python >= 3.2 podrias usar también itertools.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 usar np.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 simplemente:
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 dentro taobeta
y taogamma
porque usas itertools.starmap y no es posible pasar los argumentos b
,n
,y
y s
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 de calcular
. Nota que la segunda definición de p
que haces en tu código (la que usa taobeta
) ahora se llama p2
.
Esto son solo ideas por si te interesa aplicarlas. Seguramente se pueden hacer otras mejoras pero hay conceptos del procedimiento matemático de tu código que lógicamente se me escapan por lo que no me atrevo a tocar más sin saber el significado de cada variable/función.
for
solo actualizas B, ¿deseas que tambien se vean reflejados los cambios enLA
,y
,LB
yb
según cambia el contenido deB
?B
en el loop que lo actualice en la parte inicial del codigoLA
LB
para hacer que me itere todo el codigo nuevamente Mil gracias FJsevilla.B
las demás no se modifican porque son objetos distintos, son copias de los datos deB
pero no van a cambiar. Si deseas solo modicar esos valores añade alfor
las mismas declaraciones deLA
,y
LB
yb
para que se actualicen todas. ¿O acaso quieres que todo el código despues de declararB
,LA
,LB
,y
yb
hasta elfor
se ejecute de nuevo con el nuevo valor de B? Si ese es el caso encierralo en una función y en cada iteracion del for la llamas.B
,LA
,LB
,y
yb
hasta el for se ejecute de nuevo con el nuevoB
Lo que dices es que la parte inicial no el loop final, lo encierre en un def _() con un returny
yb
Si me puedes dar un poco mas de luz porque llevo dias con este problema, he visto otros codigos pero no logro hacer que este me corra como deberia. @FJSevilla