Skip to main content
2 de 4
se añadieron 413 caracteres en el cuerpo
FJSevilla
  • 62.1k
  • 7
  • 43
  • 68

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.

FJSevilla
  • 62.1k
  • 7
  • 43
  • 68