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. El código en general se puede hacer más eficiente, por ejemplo las funciones (p
, UI
, t
, etc) deberian estar declaradas fuera de calcular()
para no volver a declararlas en cada iteración. No quiero modificar demasiado tu código sin saber que hace cada cosa pero hay algunas cosas que no termino de entender, por ejemplo a la función UI necesita 5 parámetros de entrada pero solo usas uno dentro de ella.