Skip to main content
4 de 4
se añadieron 4697 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.

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.

FJSevilla
  • 62.1k
  • 7
  • 43
  • 68