|
|
< >一元线性回归程序，优点在于可得出拟合优度概率等参数，刚从一个java程序改过来的，不知道有没有改出错来。请大家帮忙测试一下。</P>
<>#include <math.h>
#include <stdlib.h>
#include <iostream.h></P>
<>class CDataFit
{
public:
double fit_a;//y=a+bx,存放系数 a 。
double fit_b;//y=a+bx,存放系数 b 。
double fit_siga;//a 的标准差。
double fit_sigb;//b 的标准差。
double fit_chi2;//sum[（yi-a-bxi）^2/delta^2]
double fit_q;//拟合优度概率，当mwt=0时，fit_q=1 </P>
<P>void fit(double x[], /*存放自变量的ndata个观测值，注意从x[1]开始*/
double y[], /*存放 y 的ndata个观测值。注意从y[1]开始*/
int ndata, /*数据量 */
double sig[], // y 的测量误差的标准差。
int mwt /* mwt=0，则sig[]中的数据无效，不予考虑。*/
);
private:
double gammq(double a, double x);//计算不完全gammar函数
void gser(double &gamser, double , double x, double &gln);//gammp内部调用
void gcf(double &gammcf, double a, double x, double &gln);//gammp内部调用
double gammln(double xx);//计算gammar[xx]的对数，即Ln[gammar[xx]]
};</P>
<P>void CDataFit::fit(double x[],double y[],int ndata,double sig[],int mwt)
{
double a= 0.0;
double b= 0.0;
double siga= 0.0;
double sigb= 0.0;
double chi2= 0.0;
double q= 0.0;
int i;
double sigdat, t, sxoss, wt, ss, sx = 0.0;
double sy = 0.0;
double st2 = 0.0;
if (mwt != 0)
{
ss = 0.0;
for (i = 1; i <= ndata; i++)
{
wt = 1.0 / (sig * sig);
ss = ss + wt;
sx = sx + x * wt;
sy = sy + y * wt;
}
}
else
{
for (i = 1; i <= ndata; i++)
{
sx = sx + x;
sy = sy + y;
}
ss = ndata;
}
sxoss = sx / ss;
if (mwt != 0)
{
for (i = 1; i <= ndata; i++)
{
t = (x - sxoss) / sig;
st2 = st2 + t * t;
b = b + t * y / sig;
}
}
else
{
for (i = 1; i <= ndata; i++)
{
t = x - sxoss;
st2 = st2 + t * t;
b = b + t * y;
}
}
b = b / st2;
a = (sy - sx * b) / ss;
siga = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
sigb = sqrt(1.0 / st2);
chi2 = 0.0;
if (mwt == 0)
{
for (i = 1; i <= ndata; i++)
{
chi2 = chi2 + (y - a - b * x)*(y - a - b * x);
}
q = 1.0;
sigdat = sqrt(chi2 / (ndata - 2));
siga = siga * sigdat;
sigb = sigb * sigdat;
}
else
{
for (i = 1; i <= ndata; i++)
{
chi2 = chi2 + pow(((y - a - b * x) / sig) , 2.0);
}
q = gammq(0.5 * (ndata - 2), 0.5 * chi2);
}
fit_a = a; fit_b = b; fit_siga = siga;
fit_sigb = sigb; fit_chi2 = chi2; fit_q = q;
}</P>
<P>
double CDataFit::gammq(double a, double x)
{
double temp, gamser, gammcf, gln;
gamser = 0; gln = 0; gammcf = 0;
if (x < 0.0 || a <= 0.0)
{
exit(1);
}
if (x < a + 1.0)
{
gser(gamser, a, x, gln);
temp= 1.0 - gamser;
}
else
{
gcf(gammcf, a, x, gln);
temp = gammcf;
}
return temp;
}</P>
<P>void CDataFit::gser(double &gamser, double a, double x, double &gln)
{
int itmax, n;
double ap, sum, del, eps;
itmax = 100;
eps = 0.0000003;
gln = gammln(a);
if (x <= 0.0)
{
if (x < 0.0)
{
exit(1);
}
gamser = 0.0;
return;
}
ap = a;
sum = 1.0 / a;
del = sum;
for( n = 1; n <= itmax; n++)
{
ap = ap + 1.0;
del = del * x / ap;
sum = sum + del;
if (fabs(del) < (fabs(sum) * eps))
{
gamser = sum * exp(-x + a * log(x) - gln);
break;
}
}
}
void CDataFit::gcf(double &gammcf, double a, double x, double &gln)
{
int itmax, n;
double eps, a0, a1, b0, b1, fac, an, ana, anf, gold, g;
itmax = 100;
eps = 0.0000003;
gln = gammln(a);
gold = 0.0;
a0 = 1.0;
a1 = x;
b0 = 0.0;
b1 = 1.0;
fac = 1.0;
for (n = 1; n <= itmax; n++)
{
an = n;
ana = an - a;
a0 = (a1 + a0 * ana) * fac;
b0 = (b1 + b0 * ana) * fac;
anf = an * fac;
a1 = x * a0 + anf * a1;
b1 = x * b0 + anf * b1;
if (a1 != 0.0)
{
fac = 1.0 / a1;
g = b1 * fac;
if (fabs((g - gold) / g) < eps )
{
gammcf = exp(-x + a * log(x) - gln) * g;
break;
}
gold = g;
}
}
}
double CDataFit::gammln(double xx)
{
int j;
double temp;
double cof[7] = {
0.0,
76.18009173,
-86.50532033,
24.01409822,
-1.231739516,
0.00120858003,
-0.00000536382,
};
double stp, half, one, fpf, x, tmp, ser;</P>
<P>stp = 2.50662827465;
half = 0.5;
one = 1.0;
fpf = 5.5;</P>
<P>x = xx - one;
tmp = x + fpf;
tmp = (x + half) * log(tmp) - tmp;
ser = one;
for (j = 1; j <= 6; j++)
{
x = x + one;
ser = ser + cof[j] / x;
}
temp = tmp + log(stp * ser);
return temp;
}</P>
<P>void main (void)//验证程序
{
int i, mwt, npt = 100;
double spread = 20;
double x[101];
double y[101];
double sig[101];
double a, b, siga, sigb, chi2, q;
CDataFit g;</P>
<P>a = 0; b = 0; siga = 0; sigb = 0; chi2 = 0; q = 0;
for (i = 1; i <= npt; i++)
{
x = 0.1 * i;
y = -2.0 * x + 1.0 + spread * rand() / 32767.;
sig = spread;
}
g.fit(x, y, npt, sig, mwt);
a = g.fit_a;
b = g.fit_b;
siga = g.fit_siga;
sigb = g.fit_sigb;
chi2 = g.fit_chi2;
q = g.fit_q;
cout<<"Ignoring standard deviation"<<endl;
cout<<"a = "<<a<<" "<<endl;
cout<<"Uncertainty: "<<siga<<endl;
cout<<"b = "<<b<<" "<<endl;
cout<<"Uncertainty: "<<sigb<<endl;
cout<<"Chi-squared: "<<chi2<<endl;
cout<<"Goodness-of-fit: "<<q<<endl;
}</P>
<P>
</P> |
|
发表于 2004-8-7 07:31:02