tp43.py (793B)
1 #!/usr/bin/env python3.9 2 3 from sympy import * 4 5 def main(): 6 p_, n, nprime = symbols('p_, n, nprime') 7 vr = Matrix([p_, n, nprime]) # array of variables 8 9 p = Matrix([1, 0.8, 1.2]) # initial guess 10 y = Matrix([0.8, 1.2, 1, 1]) # observables 11 P = eye(4)*25 # inverse convariace matrix 12 f = Matrix(([n*p_, nprime*p_, n, nprime])) # parameters 13 14 15 for _ in range(10): 16 17 f0 = f.subs({p_:p[0], n:p[1], nprime:p[2]}) 18 I0 = y - f0 19 20 # design matrix 21 A = f.jacobian(vr).subs({p_:p[0], n:p[1], nprime:p[2]}) 22 23 dp=(A.T @ P @ A).inv() @ A.T @ P @ I0 24 25 p = p + dp 26 27 cov = (A.T @ P @ A).inv() 28 for i in range(3): 29 print(f'{vr[i]}:\t {p[i]}\t\\pm {sqrt(cov[i, i])}') 30 31 if __name__ == "__main__": 32 main()