最小二乘法曲线拟合y=c1+c2x+c3x^2——java版
程序员文章站
2022-03-26 20:35:29
public class Nllsf { /** * 拟合目标P(x)=c1+c2x+c3x^2 * c1+c2x+c3x^2=f(x) */static double[][] ATA;static double[] ATb;static double[][] AT;//A的转置static double[][] A;static double[] b;static double[] c;static void CreatAb(int n,double[][]...
import java.util.Scanner;
public class Nllsf {
/**
* 拟合目标P(x)=c1+c2x+c3x^2
* c1+c2x+c3x^2=f(x)
*/
static double[][] ATA;
static double[] ATb;
static double[][] AT;//A的转置
static double[][] A;
static double[] b;
static double[] c;
static void CreatAb(int n,double[][] xy){
A=new double[n][3];
b=new double[n];
c=new double[3];
for(int i=0;i<n;i++){
A[i][0]=Math.pow(xy[0][i],0);
A[i][1]=Math.pow(xy[0][i],1);
A[i][2]=Math.pow(xy[0][i],2);
b[i]=xy[1][i];
}
}
static void CreatATAandATbandAT(int n){
AT=new double[3][n];
for(int i=0;i<3;i++){
for(int j=0;j<n;j++){
AT[i][j]=A[j][i];
}
}
ATA=new double[3][3];
for(int i=0;i<3;i++){
for(int ii=0;ii<3;ii++){
double sum=0;
for(int j=0;j<n;j++){
sum=AT[i][j]*AT[ii][j]+sum;
ATA[i][ii]=sum;
}
}
}
ATb=new double[3];
for(int i=0;i<3;i++){
double sum=0;
for(int j=0;j<n;j++){
sum=sum+AT[i][j]*b[j];
}
ATb[i]=sum;
}
}
static void Guss(){ //高斯列主元消去法
double[][] ATAATb=new double[3][4];
int num=ATAATb[0].length-1;
for(int i=0;i<num;i++){
for(int j=0;j<num+1;j++){
if(j==num){
ATAATb[i][j]=ATb[i];
}
else{
ATAATb[i][j]=ATA[i][j];
}
}
}
for(int ij=0;ij<num-1;ij++){//高斯类主元消去法突破点:主元素
double max=ATAATb[ij][ij];
int maxi=ij;
for(int i=ij+1;i<num;i++){//找到主元素所在列的最大值
if(ATAATb[i][ij]>=max){
max=ATAATb[i][ij];
maxi=i;
}
}
for(int i=ij+1;i<num;i++){//清零
double k=ATAATb[i][ij]/ATAATb[ij][ij];
for(int j=ij;j<num+1;j++){
ATAATb[i][j]=ATAATb[i][j]-ATAATb[ij][j]*k;
}
}
}
//以上为高斯列主元消去法的准备工作
c[num-1]=ATAATb[num-1][num]/ATAATb[num-1][num-1];
for(int i=num-2;i>-1;i--){
double sum=ATAATb[i][num];
for(int j=num-1;j>i;j--){
sum=sum-ATAATb[i][j]*c[j];
}
c[i]=sum/ATAATb[i][i];
}//回代求值
}
public static void main(String[] args) {
Scanner in=new Scanner(System.in);
int n=in.nextInt();
double[][] xy=new double[2][n];//xy[0][]放x,xy[1][]放y
for(int i=0;i<n;i++){
xy[0][i]=in.nextDouble();
xy[1][i]=in.nextDouble();
}
CreatAb(n, xy);
CreatATAandATbandAT(n);
Guss();//高斯列主元消去法解方程组
System.out.print("Y=");
System.out.printf("%.6f",c[0]);
System.out.print("+");
System.out.printf("%.6f",c[1]);
System.out.print("X+");
System.out.printf("%.6f",c[2]);
System.out.print("X^2");
}
}
例子:
x | 1 | 3 | 4 | 5 | 6 | 7 | 8 |
---|---|---|---|---|---|---|---|
y | 10 | 5 | 4 | 2 | 1 | 1 | 2 |
本文地址:https://blog.csdn.net/Ken17963/article/details/107598425