#------------------- polylsq.g ---------------------------- # # Glish script to demonstrate using matrix object to do # least squares fitting to a polynomial with error analysis # #tablefromascii('begtable','begdata.txt','begdatahdr') xydata := table('begtable') xydata.getkeywords() # Print out needed keywords x := xydata.getcol('XDATA') y := xydata.getcol('YDATA') N := length(y) ydata := array(y,N,1) # Need array not vector for matrix input A := array(0,N,3) A[,1] := 1 for (i in 1:N) for (j in 2:3) A[i,j] := x[i]^(j-1) mx := matrix_functions() Atran := mx.transpose(A) G := mx.mult(Atran,A) Ginv := mx.invert(G) Coef := mx.mult(Ginv,mx.mult(Atran,ydata)) # Solve for polynomial coefficients ysoln := mx.mult(A,Coef) # Polynomial fit to y(x) errsq := N*mx.mean((ydata - ysoln)^2)/(N-1) Cov := errsq*Ginv # Covariance matrix rms := errsq^0.5 for (j in 1:3) sigma[j] := Cov[j,j]^0.5 # Errors for polynomial coefficients r := Cov # Compute r = Correlation Matrix for (j in 1:3) for (k in 1:3) r[j,k] := Cov[j,k]/(Cov[j,j]*Cov[k,k])^0.5 pl := pgplotter(background="white") # Set up plot to plot data and fit red := 2 blue := 4 pl.sci(blue) pl.env(min(x),max(x),min(y),max(y),0,1) pl.lab("x","y","Plot Data (Red) and 2nd Order Polynomial Fit to Data (Blue)") pl.line(x,ysoln) pl.sci(red) pl.pt(x,ydata,4)