0

Tengo el siguiente codigo, necesito que b, y sean actualizados intente un for in al final del codigo para el actualizara los valores del inicio pero no dio resultado un buen resultado.

Code

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]])

Esta es la parte del codigo que se deberia actualizar con el loop al final

B = np.matrix([[0.5673],[0.0026]])
LA= np.array(np.array(B)[:,0])
y= np.array(LA[0])
LB= np.array(B[1])[0]
b= np.array(LB[0])

La parte del codigo que lo actualiza esta al final

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]]])

Esta es la parte del codigo el for in queria que le enviara el dato a la parte inicial del codigo

i=1
for _ in range(2):

   i=i+1

   B= (dot(inv(jacobiana.T*jacobiana),jacobiana.T*(iteorica-itecho)))+B

El problema es que no actualiza los parametros superiores.

5
  • En el for solo actualizas B, ¿deseas que tambien se vean reflejados los cambios en LA, y, LB y b según cambia el contenido de B?
    – FJSevilla
    Commented el 25 mar. 2017 a las 15:58
  • Si asi es, al cambiar el contenido de B en el loop que lo actualice en la parte inicial del codigo LA LB para hacer que me itere todo el codigo nuevamente Mil gracias FJsevilla. Commented el 25 mar. 2017 a las 16:09
  • Ten en cuenta que los enteros o floats son objetos inmutables, si modificas B las demás no se modifican porque son objetos distintos, son copias de los datos de B pero no van a cambiar. Si deseas solo modicar esos valores añade al for las mismas declaraciones de LA, y LB y b para que se actualicen todas. ¿O acaso quieres que todo el código despues de declarar B, LA, LB, y y b hasta el for 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.
    – FJSevilla
    Commented el 25 mar. 2017 a las 16:27
  • Si eso es lo intento que todo el codigo despues de declarar B, LA, LB, y y b hasta el for se ejecute de nuevo con el nuevo B Lo que dices es que la parte inicial no el loop final, lo encierre en un def _() con un return y y b 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 Commented el 25 mar. 2017 a las 16:42
  • Vale, creo que te he entendido... Voy a crear una respuesta a ver si es lo que deseas, sino vamos trabajando sobre ella.
    – FJSevilla
    Commented el 25 mar. 2017 a las 16:43

1 respuesta 1

2

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.

7
  • Primero decirte que apreio mucho esta ayuda. El codigo corre, pero tengo un poblema con print('Iteración {}, valor de B:\n {}'.format(i, B)) ya que con esa linea de comando el codigo no corre, estoy usando python 27.1 y no se si tenga algo que ver o sea alguna libreria que necesito importar? Por lo que me dices de la funcion UI como seria con los 5 parametros? Lo de numpy ya lo voy a quitar. Tienes mas observaciones del codigo te lo agradezco. Commented el 25 mar. 2017 a las 19:28
  • ¿Que ocurre, te marca algún error? Lo único que puede pasarte es que no codifique la tilde si estas en Python 2.x , de todas formas eliminala si quieres, el print es solo informativo para mostrar el valor de B en cada iteracion. O prueba con print('Iteración {}, valor de B:\n {}'.format(i, B).decode('utf8'))
    – FJSevilla
    Commented el 25 mar. 2017 a las 19:48
  • Efectivamente, era la tilde. Sin querer abusar de tu ayuda, Como seria la funcion de UI con los 5 parametros? Commented el 25 mar. 2017 a las 19:56
  • @SantiagoYepes no te preocupes por abusar, estamos para ayudar y aprender entre todos. He echado un vistaso a tu código y he editado la respuesta con algunas ideas de mejora por si te interesa. Saludos.
    – FJSevilla
    Commented el 25 mar. 2017 a las 23:17
  • Wo, FJS, muchisimas gracias!!!! voy a revisar bien lo que pones para aplicarlas. al codigo, Muchisimas gracias. Commented el 26 mar. 2017 a las 14:48

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.