📄 d9r10f.java
字号:
public class d9r10F
{
double aa,abdevt;
int ndatat;
double xt[] = new double[1001];
double yt[] = new double[1001];
double arr[] = new double[1001];
double medift_a;
double medift_b;
double medift_abdev;
void medfit(double x[], double y[], int ndata, double a, double b, double abdev)
{
int j;
double bb,del,chisq,f,sx;
double sigb,b1,b2,f1,f2,sy = 0.0;
double sxy = 0.0;
double sxx = 0.0;
sx=0;
for (j = 1; j <= ndata; j++)
{
xt[j] = x[j];
yt[j] = y[j];
sx = sx + x[j];
sy = sy + y[j];
sxy = sxy + x[j] * y[j];
sxx = sxx + x[j] * x[j];
}
ndatat = ndata;
del = ndata * sxx - sx * sx;
aa = (sxx * sy - sx * sxy) / del;
bb = (ndata * sxy - sx * sy) / del;
chisq = 0.0;
for (j = 1; j <= ndata; j++)
{
chisq = chisq + Math.pow((y[j] - (aa + bb * x[j])) , 2);
}
sigb = Math.sqrt(chisq / del);
b1 = bb;
f1 = rofunc(b1);
b2 = bb + Math.abs(3 * sigb) * sgn(f1);
f2 = rofunc(b2);
while (f1 * f2 > 0.0)
{
bb = 2.0 * b2 - b1;
b1 = b2;
f1 = f2;
b2 = bb;
f2 = rofunc(b2);
}
sigb = 0.01 * sigb;
while (Math.abs(b2 - b1) > sigb)
{
bb = 0.5 * (b1 + b2);
if (bb == b1 || bb == b2)
{
break;
}
f = rofunc(bb);
if (f * f1 >= 0.0)
{
f1 = f;
b1 = bb;
}
else
{
f2 = f;
b2 = bb;
}
}
a = aa;
b = bb;
abdev = abdevt / ndata;
medift_a = a; medift_b = b; medift_abdev = abdev;
}
double rofunc(double b)
{
int n1, j, ndata;
double nml, nmh, d;
n1 = ndatat + 1;
nml = (double)n1 / 2;
nmh = n1 - nml;
for (j = 1; j <= ndatat; j++)
{
arr[j] = yt[j] - b * xt[j];
}
sort(ndatat, arr);
aa = 0.5 * (arr[(int)nml] + arr[(int)nmh]);
double sum = 0.0;
abdevt = 0.0;
for (j = 1; j <= ndatat; j++)
{
d = yt[j] - (b * xt[j] + aa);
abdevt = abdevt + Math.abs(d);
sum = sum + xt[j] * sgn(d);
}
return sum;
}
static int iset;
static double gset;
int gasdev_idum;
double gasdev(int idum)
{
double t, v1, v2, fac, r;
if (iset == 0.0)
{
do
{
v1 = 2.0 * ran1(idum) - 1.0;
idum = ran1_idum;
v2 = 2.0 * ran1(idum) - 1.0;
idum = ran1_idum;
r = v1 * v1 + v2 * v2;
}while ((r >= 1.0) || (r == 0.0));
fac = Math.sqrt(-2.0 * Math.log(r) / r);
gset = v1 * fac;
t = v2 * fac;
iset = 1;
gasdev_idum = idum;
return t;
}
else
{
t = gset;
gasdev_idum = idum;
iset = 0;
return t;
}
}
int ran1_idum;
static double r[] = new double[98];
static int ix1, ix2, ix3;
static int iff;
double ran1(int idum)
{
double rm1, rm2 ,t;
int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;
m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
m3 = 243000; ia3 = 4561; ic3 = 51349;
if ((idum < 0) || (iff == 0))
{
iff = 1;
ix1 = (ic1 - idum) % m1;
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = ix1 % m2;
ix1 = (ia1 * ix1 + ic1) % m1;
ix3 = ix1 % m3;
for (j = 1; j <= 97; j++)
{
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
}
idum = 1;
}
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
ix3 = (ia3 * ix3 + ic3) % m3;
j = 1 + (int)((97 * ix3) / m3);
if ((j > 97) || (j < 1))
{
System.out.println( "abnormal exit");
System.exit(1);
}
t = r[j];
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
ran1_idum = idum;
return t;
}
double fit_a;
double fit_b;
double fit_siga;
double fit_sigb;
double fit_chi2;
double fit_q;
void fit(double x[], double y[], int ndata, double sig[], int mwt, double a,
double b, double siga, double sigb, double chi2, double q)
{
int i;
double sigdat, t, sxoss, wt, ss, sx = 0.0;
double sy = 0.0;
double st2 = 0.0;
b = 0.0;
if (mwt != 0)
{
ss = 0.0;
for (i = 1; i <= ndata; i++)
{
wt = 1.0 / (sig[i] * sig[i]);
ss = ss + wt;
sx = sx + x[i] * wt;
sy = sy + y[i] * wt;
}
}
else
{
for (i = 1; i <= ndata; i++)
{
sx = sx + x[i];
sy = sy + y[i];
}
ss = ndata;
}
sxoss = sx / ss;
if (mwt != 0)
{
for (i = 1; i <= ndata; i++)
{
t = (x[i] - sxoss) / sig[i];
st2 = st2 + t * t;
b = b + t * y[i] / sig[i];
}
}
else
{
for (i = 1; i <= ndata; i++)
{
t = x[i] - sxoss;
st2 = st2 + t * t;
b = b + t * y[i];
}
}
b = b / st2;
a = (sy - sx * b) / ss;
siga = Math.sqrt((1.0 + sx * sx / (ss * st2)) / ss);
sigb = Math.sqrt(1.0 / st2);
chi2 = 0.0;
if (mwt == 0)
{
for (i = 1; i <= ndata; i++)
{
chi2 = chi2 + Math.pow((y[i] - a - b * x[i]), 2);
}
q = 1.0;
sigdat = Math.sqrt(chi2 / (ndata - 2));
siga = siga * sigdat;
sigb = sigb * sigdat;
}
else
{
for (i = 1; i <= ndata; i++)
{
chi2 = chi2 + Math.pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
}
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;
}
double gammq_a, gammq_x;
double gammq(double a, double x)
{
double temp, gamser, gammcf, gln;
gamser = 0; gln = 0; gammcf = 0;
if (x < 0.0 || a <= 0.0)
{
System.out.println("pause in gammq");
System.exit(1);
}
if (x < a + 1.0)
{
gser(gamser, a, x, gln);
gamser = gser_gamser; a = gser_a; x = gser_x; gln = gser_gln;
temp= 1.0 - gamser;
}
else
{
gcf(gammcf, a, x, gln);
gammcf = gcf_gammcf; a = gcf_a; x = gcf_x; gln = gcf_gln;
temp = gammcf;
}
gammq_a = a;
gammq_x = x;
return temp;
}
double gcf_gammcf, gcf_a, gcf_x, gcf_gln;
void 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 (Math.abs((g - gold) / g) < eps )
{
gammcf = Math.exp(-x + a * Math.log(x) - gln) * g;
break;
}
gold = g;
}
}
gcf_gammcf = gammcf; gcf_a = a; gcf_x = x; gcf_gln = gln;
}
double gser_gamser, gser_a, gser_x, gser_gln;
void 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)
{
System.out.println("pause in gser");
System.exit(1);
}
gamser = 0.0;
gser_gamser = gamser; gser_a = a; gser_x = x;gser_gln = gln;
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 (Math.abs(del) < (Math.abs(sum) * eps))
{
gamser = sum * Math.exp(-x + a * Math.log(x) - gln);
break;
}
}
gser_gamser = gamser; gser_a = a; gser_x = x; gser_gln = gln;
}
double gammln(double xx)
{
int j;
double temp;
double cof[] = new double[7];
double stp, half, one, fpf, x, tmp, ser;
cof[1] = 76.18009173;
cof[2] = -86.50532033;
cof[3] = 24.01409822;
cof[4] = -1.231739516;
cof[5] = 0.00120858003;
cof[6] = -0.00000536382;
stp = 2.50662827465;
half = 0.5;
one = 1.0;
fpf = 5.5;
x = xx - one;
tmp = x + fpf;
tmp = (x + half) * Math.log(tmp) - tmp;
ser = one;
for (j = 1; j <= 6; j++)
{
x = x + one;
ser = ser + cof[j] / x;
}
temp = tmp + Math.log(stp * ser);
return temp;
}
void sort(int n, double ra[])
{
int i, j, l, ir;
double rra;
l = (int)(n / 2) + 1;
ir = n;
do
{
if (l > 1)
{
l = l - 1;
rra = ra[l];
}
else
{
rra = ra[ir];
ra[ir] = ra[1];
ir = ir - 1;
if (ir == 1)
{
ra[1] = rra;
return;
}
}
i = l;
j = l + l;
while (j <= ir)
{
if (j < ir)
{
if (ra[j] < ra[j + 1])
{
j += 1;
}
}
if (rra < ra[j])
{
ra[i] = ra[j];
i = j;
j += j;
}
else
{
j = ir + 1;
}
}
ra[i] = rra;
}while(true);
}
int sgn(double pa)
{
if (pa > 0.0)
{
return 1;
}
else
{
if (pa < 0.0)
{
return -1;
}
}
return 0;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -