# Demo entry 6645807

py

Submitted by 张樱潇 on Oct 11, 2017 at 17:50
Language: Fortran. Code size: 4.1 kB.

```program zyx04
implicit none
real,dimension(10)::x1,x2,x3,x4,y,x0,y0    !x0,y0是原始自变量，x1,x2,x3,x4,y为变换后的自变量，转换为多元线性拟合
integer,parameter::n=4                     !n为变量的个数
real t,t1,t2,t3,t4                         !分别为各项代换后变量的平均值
real b0                                    !常系数
real::s(n,n)=0.0                           !xi,xj之间的协方差矩阵
real,dimension(n)::ss=0.0,b                !ss为xi,y之间的协方差矩阵，b为系数矩阵
integer i,j
real::w=0.0                                   !求出拟合函数与原始点之间的误差

x0(1:10)=(/1.02,0.95,0.87,0.77,0.67,0.56,0.44,0.30,0.16,0.01/)
y0(1:10)=(/0.39,0.32,0.27,0.22,0.18,0.15,0.13,0.12,0.13,0.15/)

do i=1,10                                  !为每个新定义变量的数组赋值
x1(i)=x0(i)
x2(i)=y0(i)
x3(i)=x0(i)*y0(i)
x4(i)=y0(i)*y0(i)
y(i)=x0(i)*x0(i)
end do

t=0
t1=0
t2=0
t3=0
t4=0

do i=1,10                                   !求新定义变量的平均值
t1=t1+x1(i)
t2=t2+x2(i)
t3=t3+x3(i)
t4=t4+x4(i)
t=t+y(i)
end do

t1=t1/10
t2=t2/10
t3=t3/10
t4=t4/10
t=t/10

do j=1,10                                   !求xi,xj之间的协方差
s(1,1)=s(1,1)+(x1(j)-t1)*(x1(j)-t1)
s(1,2)=s(1,2)+(x1(j)-t1)*(x2(j)-t2)
s(1,3)=s(1,3)+(x1(j)-t1)*(x3(j)-t3)
s(1,4)=s(1,4)+(x1(j)-t1)*(x4(j)-t4)
s(2,2)=s(2,2)+(x2(j)-t2)*(x2(j)-t2)
s(2,3)=s(2,3)+(x2(j)-t2)*(x3(j)-t3)
s(2,4)=s(2,4)+(x2(j)-t2)*(x4(j)-t4)
s(3,3)=s(3,3)+(x3(j)-t3)*(x3(j)-t3)
s(3,4)=s(3,4)+(x3(j)-t3)*(x4(j)-t4)
s(4,4)=s(4,4)+(x4(j)-t4)*(x4(j)-t4)
end do
s(2,1)=s(1,2)
s(3,1)=s(1,3)
s(4,1)=s(1,4)
s(3,2)=s(2,3)
s(4,2)=s(2,4)
s(4,3)=s(3,4)

do i=1,4
do j=1,4
s(i,j)=s(i,j)/10
end do
end do

do i=1,10                                  !求xi,y之间的协方差
ss(1)=ss(1)+(x1(i)-t1)*(y(i)-t)
ss(2)=ss(2)+(x2(i)-t2)*(y(i)-t)
ss(3)=ss(3)+(x3(i)-t3)*(y(i)-t)
ss(4)=ss(4)+(x4(i)-t4)*(y(i)-t)
end do

do i=1,4
ss(i)=ss(i)/10
end do

write(*,*)"协方差矩阵"
do i=1,4
write(*,*)s(i,:)
end do

write(*,*)"协方差逆矩阵"
call ni(s,n)                              !调用函数，求s的逆矩阵
do i=1,4
write(*,*)s(i,:)
end do

b=matmul(s,ss)
b0=t-t1*b(1)-t2*b(2)-t3*b(3)-t4*b(4)

write(*,*)"系数"
write(*,*)b0,b

do i=1,10
w=w+(b0+b(1)*x1(i)+b(2)*x2(i)+b(3)*x3(i)+b(4)*x4(i)-y(i))**2
end do

write(*,*)"误差"
write(*,*)w

open(1,file='4.txt',form='formatted',status='replace')
write(1,*)"系数"
write(1,*)b0,b
write(1,*)"误差"
write(1,*)w
close(1)

stop
end program

subroutine ni(A,N)    !求解逆矩阵的子程序，子程序参考徐良士算法集brinv逆矩阵求解算法
implicit none
integer::N
integer::I,J,K
real::A(N,N),IS(N),JS(N)
real::T,D,L
L=1
DO 100 K=1,N
D=0.0
DO 10 I=K,N
DO 10 J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
IS(K)=I
JS(K)=J
END IF
10	  CONTINUE
IF (D+1.0.EQ.1.0) THEN
L=0
WRITE(*,20)
RETURN
END IF
20	  FORMAT(1X,'ERR**NOT INV')
DO 30 J=1,N
T=A(K,J)
A(K,J)=A(IS(K),J)
A(IS(K),J)=T
30	  CONTINUE
DO 40 I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
40	  CONTINUE
A(K,K)=1/A(K,K)
DO 50 J=1,N
IF (J.NE.K) THEN
A(K,J)=A(K,J)*A(K,K)
END IF
50	  CONTINUE
DO 70 I=1,N
IF (I.NE.K) THEN
DO 60 J=1,N
IF (J.NE.K) THEN
A(I,J)=A(I,J)-A(I,K)*A(K,J)
END IF
60	      CONTINUE
END IF
70	  CONTINUE
DO 80 I=1,N
IF (I.NE.K) THEN
A(I,K)=-A(I,K)*A(K,K)
END IF
80	  CONTINUE
100	CONTINUE
DO 130 K=N,1,-1
DO 110 J=1,N
T=A(K,J)
A(K,J)=A(JS(K),J)
A(JS(K),J)=T
110	  CONTINUE
DO 120 I=1,N
T=A(I,K)
A(I,K)=A(I,IS(K))
A(I,IS(K))=T
120	  CONTINUE
130	CONTINUE
RETURN
END
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.