Update parameters b, and with loop

If I have understood the problem you wish, giving an initial value of B , recalculate it x times. To calculate B you use your own B , jacobiana , itecho e iteorica . What you can do is enclose the code in a function that has as input parameter B and return iteorica and jacobiana, calculate B and then call it again with the new value of B , something like this:

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

It is also possible to do it recursively.

EDITING:

As I mentioned there are some things that can be improved, here I leave the code with some changes that I will comment later:

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

The most significant changes are:

• Clean the import a bit as I mentioned, it is especially important to avoid from module import * , especially with complex and extensive modules such as NumPy or SciPy . You can end up rewriting functions of these modules unintentionally and with unpleasant results in general.

• In many of your operations it is plausible to use list-comprehensions avoiding unnecessary variables and improving efficiency, for example:

  fi=[]
  for x in Xu:
      qq=ii(x, b, n, y, s)
      fi.append(qq)
  

  It may simply be like:

  fi=[ii(x, b, n, y, s) for x in Xu]
  

• You make several accumulated sums using cycles for , you can do this efficiently using numpy.cumsum If you were using a version of Python > = 3.2 you could also use itertools.acumulate . In this way you can efficiently simplify things like:

  suma=0
  sumtaobeta=[]
  for i in valtaobeta:
      suma=suma+i
      sumtaobeta.append(suma)
  

  By:

  sumtaobeta = np.cumsum(valtaobeta)
  

• To calculate jacobiana you can use np.column_stack , reducing:

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

  to simply:

  jacobiana = np.matrix(np.column_stack((multib,g)))
  

• All the functions that are easy to get out of calcular() are now out, only taobeta and taogamma remain because you use itertools.starmap and it is not possible to pass the arguments b , n , y and s through it. You could also get some changes but I do not think it's worth it. This avoids that every time we call calculate to redefine all these functions, now they are global and are not within calcular . Note that the second definition of p that you make in your code (the one that uses taobeta ) is now called p2 .

These are just ideas if you are interested in applying them. Surely you can make other improvements but there are concepts of the mathematical procedure of your code that logically escape me so I dare not touch more without knowing the meaning of each variable / function.