kikiririさんへのお返事です。
Full BASICの場合、行列が計算できるので、他の言語よりは簡単に計算できると思います。
ただ、２項演算までですから、展開しながらこつこつ計算する必要があります。
また、表計算の方がGUIを含めて実用化し易いかもしれません。
上記サイトの例題（PDFファイル内）のサンプルコーディング
!最小２乗法による測量網平均（条件方程式法）
!Ｈ型
! Ａ Ｃ
! 1↓ 5 ↓3
! Ｐ → Ｑ
! 2↑ ↑4
! Ｂ Ｄ
!既知点
LET C=4 !数
DATA 25.645 !Ａ点の標高（m）
DATA 24.666 !Ｂ点
DATA 25.024 !Ｃ点
DATA 25.699 !Ｄ点
DIM Z(C)
MAT READ Z
!観測値
LET P=5 !数
DATA -6.225, 0.44 !路線１ 高低差（m）、路線長（km）
DATA -5.245, 0.25 !路線２
DATA 0.278, 0.33 !路線３
DATA -0.399, 0.26 !路線４
DATA 5.879, 0.44 !路線５
DIM X(P),G(P,P) !観測値、コアファクタ
MAT G=ZER
FOR i=1 TO P
READ X(i),G(i,i)
NEXT i
MAT PRINT X;
MAT PRINT G;
!未知数
LET Px=2 !求点ＰとＱ
!------------------------------
!自由度R ※条件方程式の数
LET R=P-Px
!条件方程式 UV=t
DIM U(R,P)
DATA 1,-1,0,0,0 !点Ｐについて HA+h1~=HB+h2~ より、(h1+v1)-(h2+v2)=HB-HA ∴v1-v2=-(h1+HA)+(h2+HB)
DATA 0,0,1,-1,0 !点Ｑについて HC+h3~=HD+h4~
DATA 1,0,0,-1,1 !点Ｐと点Ｑについて HA+h1~+h5~=HD+h4~
MAT READ U
DIM t(R,1)
FOR i=1 TO R
LET s=0
FOR j=1 TO P
IF j>C THEN LET Zj=0 ELSE LET Zj=Z(j)
LET s=s-(X(j)+Zj)*U(i,j)
NEXT j
LET t(i,1)=s
NEXT i
!!!LET t(1,1)=-X(1)+X(2)-Z(1)+Z(2) !0.001
!!!LET t(2,1)=-X(3)+X(4)-Z(3)+Z(4) !-0.002
!!!LET t(3,1)=-X(1)+X(4)-X(5)-Z(1)+Z(4) !0.001
MAT PRINT t;
!------------------------------
!相関式 NK=t、N=UG(Ut)より、K=(Ni)tを求める
DIM Ut(P,R)
MAT Ut=TRN(U)
DIM TMP(P,R) !G(Ut) ※次でも使う
MAT TMP=G*Ut
DIM N(R,R)
MAT N=U*TMP
MAT PRINT N;
DIM invN(R,R)
MAT invN=INV(N)
DIM K(R,1)
MAT K=invN*t
MAT PRINT K;
!補正値の計算（mm） V=G(Ut)K
DIM V(P,1)
MAT V=TMP*K
MAT PRINT V;
!最確値 X~=X+V
FOR i=1 TO P
PRINT "路線";STR$(i);"=";X(i)+V(i,1)
NEXT i
PRINT "点Ｐ";Z(1)+(X(1)+V(1,1)); Z(2)+(X(2)+V(2,1)) !有効桁数 ##.###
PRINT "点Ｑ";Z(3)+(X(3)+V(3,1)); Z(4)+(X(4)+V(4,1))
!精度の計算（mm^2） σ^2=(Kt)NK/r
DIM s2(1,1),Kt(1,R)
MAT Kt=TRN(K)
MAT s2=Kt*t !NK=t
LET sigma2=s2(1,1)/R
PRINT "σ^2=";sigma2
!分散行列（mm^2） σ^2*Gx、Gx=G-G(Ut)(Ni)UG
DIM Gx(P,P)
MAT Gx=TMP*invN !G(Ut)(Ni)
MAT Gx=Gx*U
MAT Gx=Gx*G
MAT Gx=G-Gx
MAT Gx=sigma2*Gx
MAT PRINT Gx;
END