FUNC1.CPP
上传用户:zengbais
上传日期:2022-08-08
资源大小:49k
文件大小:3k
开发平台:

C++ Builder

  1. #include <values.h>
  2. #include <math.h>
  3. #include <stdio.h>
  4. #include "func.h"
  5. #include "cmatrix.h"
  6. algopair::algopair(matrix& xy, size_t m, DOUBLE & dt0,DOUBLE &dt1, DOUBLE &dt2)
  7. // 最小二乘拟合构造函数
  8. // m为拟合多项式的项数,dt0为误差平方和,dt1为误差绝对值和,
  9. // dt2为最大误差绝对值
  10. {
  11. size_t n = xy.rownum;
  12. size_t i,j,k;
  13. DOUBLE z,p,c,g,q,d1,d2,s[20],t[20],b[20];
  14. if (m>n) m=n;
  15. if (m>20) m=20;
  16. data = matrix(m);
  17. data = 0.0;
  18. z=0.0;
  19. for (i=0; i<n; i++) z+=xy(i)/(1.0*n);
  20. b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
  21. for (i=0; i<=n-1; i++)
  22. {p+=xy(i)-z; c+=xy(i,1);}
  23. c/=d1; p/=d1;
  24. data.set(0,c*b[0]);
  25. if (m>1)
  26. { t[1]=1.0; t[0]=-p;
  27.   d2=0.0; c=0.0; g=0.0;
  28.   for (i=0; i<=n-1; i++)
  29.  { q=xy(i)-z-p; d2=d2+q*q;
  30. c=c+xy(i,1)*q;
  31. g=g+(xy(i)-z)*q*q;
  32.  }
  33.   c=c/d2; p=g/d2; q=d2/d1;
  34.   d1=d2;
  35.   data.set(1,c*t[1]);
  36.   data.set(0,c*t[0]+data(0));
  37. }
  38. for (j=2; j<=m-1; j++)
  39. { s[j]=t[j-1];
  40.   s[j-1]=-p*t[j-1]+t[j-2];
  41.   if (j>=3)
  42.  for (k=j-2; k>=1; k--)
  43. s[k]=-p*t[k]+t[k-1]-q*b[k];
  44.   s[0]=-p*t[0]-q*b[0];
  45.   d2=0.0; c=0.0; g=0.0;
  46.   for (i=0; i<=n-1; i++)
  47.  { q=s[j];
  48. for (k=j; k>0; k--)
  49.   q=q*(xy(i)-z)+s[k-1];
  50. d2+=q*q; c+=xy(i,1)*q;
  51. g+=(xy(i)-z)*q*q;
  52.  }
  53.   c=c/d2; p=g/d2; q=d2/d1;
  54.   d1=d2;
  55.   data.set(j,c*s[j]);
  56.   t[j]=s[j];
  57.   for (k=j; k>0; k--)
  58.  { data.set(k-1,data(k-1)+c*s[k-1]);
  59. b[k-1]=t[k-1]; t[k-1]=s[k-1];
  60.  }
  61. }
  62.  dt0=0.0; dt1=0.0; dt2=0.0;
  63.  for (i=0; i<=n-1; i++)
  64. { q=data(m-1);
  65.   for (k=m-1; k>0; k--)
  66.  q=data(k-1)+q*(xy(i)-z);
  67.   p=q-xy(i,1);
  68.   if (fabs(p)>dt2) dt2=fabs(p);
  69.   dt0+=p*p;
  70.   dt1+=fabs(p);
  71. }
  72.  xshift = 0.0;
  73.  for (i=0; i<n; i++)
  74. xshift += xy(i);
  75.  xshift /= n;
  76. }
  77. funcpair::funcpair(matrix& xy, size_t m, DOUBLE & t0,DOUBLE &t1,DOUBLE &t2):
  78. // xy为n行2列数组,存放n个样本点
  79. // m为拟合多项式的项数,dt0为误差平方和,dt1为误差绝对值和,
  80. // dt2为最大误差绝对值
  81. func(new algopair(xy,m,t0,t1,t2)) {}
  82. cmatrix algopoly::getroots() // 求出此多项式的所有根
  83. {
  84. size_t n = data.rownum-1;
  85. matrix m(n,n);
  86. matrix a,b;
  87. m = 0.0; // 先初始化成零矩阵
  88. size_t i;
  89. DOUBLE s;
  90. s = data(n);
  91. for(i=0;i<n;i++)
  92. m.set(0,n-i-1,-data(i)/s);
  93. for(i=1;i<n;i++)
  94. m.set(i,i-1,1.0);
  95. m.qreigen(a,b);
  96. return cmatrix(a,b);
  97. }
  98. algoregress::algoregress(matrix & xy, matrix* dt)
  99. // xy为n行2列的矩阵为n个样本点
  100. // dt必须指向一个6列一行的矩阵向量,六个数依次为偏差平方和,平均标准
  101. // 偏差,回归平方和,最大偏差,最小偏差,偏差平均值
  102. {
  103. size_t i,n;
  104. double xx,yy,e,f,q,u,p,umax,umin,s,a,b;
  105. xx=0.0; yy=0.0;
  106. n = xy.rownum;
  107. for (i=0; i<n; i++)
  108. { xx+=xy(i); yy+=xy(i,1);}
  109. xx/=n; yy/=n;
  110. e=0.0; f=0.0;
  111. for (i=0; i<n; i++)
  112. { q=xy(i)-xx; e+=q*q;
  113.   f+=q*(xy(i,1)-yy);
  114. }
  115. a=f/e; b=yy-a*xx;
  116. q=0.0; u=0.0; p=0.0;
  117. umax=0.0; umin=MAXDOUBLE;
  118. for (i=0; i<n; i++)
  119. { s=a*xy(i)+b;
  120.   q=q+(xy(i,1)-s)*(xy(i,1)-s);
  121.   p=p+(s-yy)*(s-yy);
  122.   e=fabs(xy(i,1)-s);
  123.   if (e>umax) umax=e;
  124.   if (e<umin) umin=e;
  125.   u+=e/n;
  126. }
  127. yfactor = a;
  128. addconst = b;
  129. if(!dt || dt->rownum < 6) return;
  130. dt->set(0,q);
  131. dt->set(1,sqrt(q/n));
  132. dt->set(2,p);
  133. dt->set(3,umax);
  134. dt->set(4,umin);
  135. dt->set(5,u);
  136. }
  137. funcregress::funcregress(matrix & xy, matrix* dt):
  138.   // xy为n行2列的矩阵为n个样本点
  139. // dt必须指向一个6列一行的矩阵向量,六个数依次为偏差平方和,平均标准
  140. // 偏差,回归平方和,最大偏差,最小偏差,偏差平均值
  141.  func(new algoregress(xy,dt)){}