C#，数值计算——数据建模Fitexy的计算方法与源程序

news2024/11/23 19:00:21

1 文本格式

using System;

namespace Legalsoft.Truffer
{
public class Fitexy
{
private double a { get; set; }
private double b { get; set; }
private double siga { get; set; }
private double sigb { get; set; }
private double chi2 { get; set; }
private double q { get; set; }
private int ndat { get; set; }
private double[] xx { get; set; }
private double[] yy { get; set; }
private double[] sx { get; set; }
private double[] sy { get; set; }
private double[] ww { get; set; }
private double aa;
private double offs;

public Fitexy(double[] x, double[] y, double[] sigx, double[] sigy)
{
this.ndat = x.Length;
this.xx = new double[ndat];
this.yy = new double[ndat];
this.sx = new double[ndat];
this.sy = new double[ndat];
this.ww = new double[ndat];
const double POTN = 1.571000;
const double BIG = 1.0e30;
const double ACC = 1.0e-6;

Gamma gam = new Gamma();
//Brent brent = new Brent(ACC);
Chixy chixy = new Chixy(xx, yy, sx, sy, ww, ref aa, ref offs);
double varx = 0.0;
double vary = 0.0;
double[] ang = new double[7];
double[] ch = new double[7];
double dum1 = 0.0;
Moment.avevar(x, ref dum1, ref varx);
Moment.avevar(y, ref dum1, ref vary);
double scale = Math.Sqrt(varx / vary);
for (int j = 0; j < ndat; j++)
{
xx[j] = x[j];
yy[j] = y[j] * scale;
sx[j] = sigx[j];
sy[j] = sigy[j] * scale;
ww[j] = Math.Sqrt(Globals.SQR(sx[j]) + Globals.SQR(sy[j]));
}
Fitab fit = new Fitab(xx, yy, ww);
b = fit.b;
offs = ang[0] = 0.0;
ang[1] = Math.Atan(b);
ang[3] = 0.0;
ang[4] = ang[1];
ang[5] = POTN;
for (int j = 3; j < 6; j++)
{
ch[j] = chixy.funk(ang[j]);
}

Brent brent = new Brent(ACC);
//brent.func = chixy.functorMethod;
brent.bracket(ang[0], ang[1], chixy);
ang[0] = brent.ax;
ang[1] = brent.bx;
ang[2] = brent.cx;
ch[0] = brent.fa;
ch[1] = brent.fb;
ch[2] = brent.fc;
b = brent.minimize(chixy);
chi2 = chixy.funk(b);
a = aa;
q = gam.gammq(0.5 * (ndat - 2), chi2 * 0.5);
double r2 = 0.0;
for (int j = 0; j < ndat; j++)
{
r2 += ww[j];
}
r2 = 1.0 / r2;
double bmx = BIG;
double bmn = BIG;
offs = chi2 + 1.0;
for (int j = 0; j < 6; j++)
{
if (ch[j] > offs)
{
double d1 = Math.Abs(ang[j] - b);
while (d1 >= Math.PI)
{
d1 -= Math.PI;
}
double d2 = Math.PI - d1;
if (ang[j] < b)
{
Globals.SWAP(ref d1, ref d2);
}
if (d1 < bmx)
{
bmx = d1;
}
if (d2 < bmn)
{
bmn = d2;
}
}
}
if (bmx < BIG)
{
bmx = Roots.zbrent(chixy, b, b + bmx, ACC) - b;
double amx = aa - a;
bmn = Roots.zbrent(chixy, b, b - bmn, ACC) - b;
double amn = aa - a;
sigb = Math.Sqrt(0.5 * (bmx * bmx + bmn * bmn)) / (scale * Globals.SQR(Math.Cos(b)));
siga = Math.Sqrt(0.5 * (amx * amx + amn * amn) + r2) / scale;
}
else sigb = siga = BIG;
a /= scale;
b = Math.Tan(b) / scale;
}
}

public class Chixy : UniVarRealValueFun
{
private double[] xx { get; set; }
private double[] yy { get; set; }
private double[] sx { get; set; }
private double[] sy { get; set; }
private double[] ww { get; set; }
private double aa { get; set; }
private double offs { get; set; }

public Chixy(double[] xxx, double[] yyy, double[] ssx, double[] ssy, double[] www, ref double aaa, ref double ooffs)
{
this.xx = xxx;
this.yy = yyy;
this.sx = ssx;
this.sy = ssy;
this.ww = www;
this.aa = aaa;
this.offs = ooffs;
}

public double funk(double bang)
{
return get(bang);
}

public double get(double bang)
{
const double BIG = 1.0e30;
//int j;
int nn = xx.Length;
//double ans;
double avex = 0.0;
double avey = 0.0;
double sumw = 0.0;
//double b;
double b = Math.Tan(bang);
for (int j = 0; j < nn; j++)
{
ww[j] = Globals.SQR(b * sx[j]) + Globals.SQR(sy[j]);
sumw += (ww[j] = (ww[j] < 1.0 / BIG ? BIG : 1.0 / ww[j]));
avex += ww[j] * xx[j];
avey += ww[j] * yy[j];
}
avex /= sumw;
avey /= sumw;
aa = avey - b * avex;
double ans = -offs;
for (int j = 0; j < nn; j++)
{
ans += ww[j] * Globals.SQR(yy[j] - aa - b * xx[j]);
}
return ans;
}
}
}

2 代码格式

using System;

namespace Legalsoft.Truffer
{
    public class Fitexy
    {
        private double a { get; set; }
        private double b { get; set; }
        private double siga { get; set; }
        private double sigb { get; set; }
        private double chi2 { get; set; }
        private double q { get; set; }
        private int ndat { get; set; }
        private double[] xx { get; set; }
        private double[] yy { get; set; }
        private double[] sx { get; set; }
        private double[] sy { get; set; }
        private double[] ww { get; set; }
        private double aa;
        private double offs;

        public Fitexy(double[] x, double[] y, double[] sigx, double[] sigy)
        {
            this.ndat = x.Length;
            this.xx = new double[ndat];
            this.yy = new double[ndat];
            this.sx = new double[ndat];
            this.sy = new double[ndat];
            this.ww = new double[ndat];
            const double POTN = 1.571000;
            const double BIG = 1.0e30;
            const double ACC = 1.0e-6;

            Gamma gam = new Gamma();
            //Brent brent = new Brent(ACC);
            Chixy chixy = new Chixy(xx, yy, sx, sy, ww, ref aa, ref offs);
            double varx = 0.0;
            double vary = 0.0;
            double[] ang = new double[7];
            double[] ch = new double[7];
            double dum1 = 0.0;
            Moment.avevar(x, ref dum1, ref varx);
            Moment.avevar(y, ref dum1, ref vary);
            double scale = Math.Sqrt(varx / vary);
            for (int j = 0; j < ndat; j++)
            {
                xx[j] = x[j];
                yy[j] = y[j] * scale;
                sx[j] = sigx[j];
                sy[j] = sigy[j] * scale;
                ww[j] = Math.Sqrt(Globals.SQR(sx[j]) + Globals.SQR(sy[j]));
            }
            Fitab fit = new Fitab(xx, yy, ww);
            b = fit.b;
            offs = ang[0] = 0.0;
            ang[1] = Math.Atan(b);
            ang[3] = 0.0;
            ang[4] = ang[1];
            ang[5] = POTN;
            for (int j = 3; j < 6; j++)
            {
                ch[j] = chixy.funk(ang[j]);
            }

            Brent brent = new Brent(ACC);
            //brent.func = chixy.functorMethod;
            brent.bracket(ang[0], ang[1], chixy);
            ang[0] = brent.ax; 
            ang[1] = brent.bx; 
            ang[2] = brent.cx;
            ch[0] = brent.fa; 
            ch[1] = brent.fb; 
            ch[2] = brent.fc;
            b = brent.minimize(chixy);
            chi2 = chixy.funk(b);
            a = aa;
            q = gam.gammq(0.5 * (ndat - 2), chi2 * 0.5);
            double r2 = 0.0;
            for (int j = 0; j < ndat; j++)
            {
                r2 += ww[j];
            }
            r2 = 1.0 / r2;
            double bmx = BIG;
            double bmn = BIG;
            offs = chi2 + 1.0;
            for (int j = 0; j < 6; j++)
            {
                if (ch[j] > offs)
                {
                    double d1 = Math.Abs(ang[j] - b);
                    while (d1 >= Math.PI)
                    {
                        d1 -= Math.PI;
                    }
                    double d2 = Math.PI - d1;
                    if (ang[j] < b)
                    {
                        Globals.SWAP(ref d1, ref d2);
                    }
                    if (d1 < bmx)
                    {
                        bmx = d1;
                    }
                    if (d2 < bmn)
                    {
                        bmn = d2;
                    }
                }
            }
            if (bmx < BIG)
            {
                bmx = Roots.zbrent(chixy, b, b + bmx, ACC) - b;
                double amx = aa - a;
                bmn = Roots.zbrent(chixy, b, b - bmn, ACC) - b;
                double amn = aa - a;
                sigb = Math.Sqrt(0.5 * (bmx * bmx + bmn * bmn)) / (scale * Globals.SQR(Math.Cos(b)));
                siga = Math.Sqrt(0.5 * (amx * amx + amn * amn) + r2) / scale;
            }
            else sigb = siga = BIG;
            a /= scale;
            b = Math.Tan(b) / scale;
        }
    }


    public class Chixy : UniVarRealValueFun
    {
        private double[] xx { get; set; }
        private double[] yy { get; set; }
        private double[] sx { get; set; }
        private double[] sy { get; set; }
        private double[] ww { get; set; }
        private double aa { get; set; }
        private double offs { get; set; }

        public Chixy(double[] xxx, double[] yyy, double[] ssx, double[] ssy, double[] www, ref double aaa, ref double ooffs)
        {
            this.xx = xxx;
            this.yy = yyy;
            this.sx = ssx;
            this.sy = ssy;
            this.ww = www;
            this.aa = aaa;
            this.offs = ooffs;
        }

        public double funk(double bang)
        {
            return get(bang);
        }

        public double get(double bang)
        {
            const double BIG = 1.0e30;
            //int j;
            int nn = xx.Length;
            //double ans;
            double avex = 0.0;
            double avey = 0.0;
            double sumw = 0.0;
            //double b;
            double b = Math.Tan(bang);
            for (int j = 0; j < nn; j++)
            {
                ww[j] = Globals.SQR(b * sx[j]) + Globals.SQR(sy[j]);
                sumw += (ww[j] = (ww[j] < 1.0 / BIG ? BIG : 1.0 / ww[j]));
                avex += ww[j] * xx[j];
                avey += ww[j] * yy[j];
            }
            avex /= sumw;
            avey /= sumw;
            aa = avey - b * avex;
            double ans = -offs;
            for (int j = 0; j < nn; j++)
            {
                ans += ww[j] * Globals.SQR(yy[j] - aa - b * xx[j]);
            }
            return ans;
        }
    }
}

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1071693.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！