【WPF学习手记】C#实现高斯拟合的模拟退火算法

算法功能：实现高斯拟合，可根据需要修改。

• 源码

public static double PeakPicker(double[] xfit, double[] yfit)
        {
            // 定义子函数 *************************
            // 数组拷贝
            void ArrCopy(double[] data, double[] DataCopy)
            {
                for (int i = 0; i < data.Length; i++)
                {
                    DataCopy[i] = data[i];
                }
            }

            // 最大值
            double ArrMax(double[] x)
            {
                double result = x[0];
                for (int i = 1; i < x.Length; i++)
                {
                    if (result < x[i])
                    {
                        result = x[i];
                    }
                }
                return result;
            }

            /// 最小值
            double ArrMin(double[] x)
            {
                double result = x[0];
                for (int i = 1; i < x.Length; i++)
                {
                    if (result > x[i])
                    {
                        result = x[i];
                    }
                }
                return result;
            }

            // 高斯拟合目标函数
            double Obj_Gaussian(double[] p, double[] x, double[] y)
            {
                double result = 0;
                for (int i = 0; i < x.Length; i++)
                {
                    result = result + Math.Pow(y[i] - p[0] * Math.Exp(-Math.Pow(x[i] - p[1], 2) / Math.Pow(p[2], 2)), 2);
                }
                return Math.Sqrt(result / y.Length);
            }

            // 条件判断
            bool IsTrue(double[] x, double[] x_min, double[] x_max)
            {
                bool result = false;
                bool Judge1 = (x[0] > x_min[0]) && (x[0] < x_max[0]);
                bool Judge2 = (x[1] > x_min[1]) && (x[1] < x_max[1]);
                bool Judge3 = (x[2] > x_min[2]) && (x[2] < x_max[2]);
                if (Judge1 && Judge2 && Judge3)
                {
                    result = true;
                }
                return result;
            }

            // 产生 0~1 之间的随机数
            double random()
            {
                var seed = Guid.NewGuid().GetHashCode();  // 这个种子很重要，不然随机数是重复的
                Random r = new Random(seed);
                int i = r.Next(0, 100000);
                return (double)i / 100000;
            }

            // 程序主体 ************************
            // 常量值
            double T0 = 1000;
            double Alpha = 0.99;
            double IteMax = 3000;

            // 变量初始化
            double[] MassComps = new double[3];
            double[] MassCompsNew = new double[3];
            double[] MassLB = new double[3];
            double[] MassUB = new double[3];
            double[] dMass = new double[3];
            MassLB[0] = ArrMin(yfit);
            MassLB[1] = ArrMin(xfit);
            MassLB[2] = 0;
            MassUB[0] = ArrMax(yfit);
            MassUB[1] = ArrMax(xfit);
            MassUB[2] = 5;

            // 初始模型
            for (int j = 0; j < 3; j++)
            {
                MassComps[j] = random() * (MassUB[j] - MassLB[j]) + MassLB[j];
            }
            // 计算目标函数
            double Fun = Obj_Gaussian(MassComps, xfit, yfit);
            // 迭代
            for (int i = 0; i < IteMax; i++)
            {
                double Tk = T0 * Math.Pow(Alpha, i);
                // 预更新
                double[] rand = new double[3];
                for (int j = 0; j < 3; j++)
                {
                    rand[j] = random();
                    dMass[j] = Tk * Math.Sign(random() - 0.5) * (Math.Pow(1 + 1.0 / Tk , Math.Abs(2 * random() - 1)) - 1);
                    MassCompsNew[j] = MassComps[j] + dMass[j] * (MassUB[j] - MassLB[j]);
                }
                // 判断
                if (IsTrue(MassCompsNew, MassLB, MassUB))
                {
                    double NewFun = Obj_Gaussian(MassCompsNew, xfit, yfit);
                    if (NewFun <= Fun)
                    {
                        ArrCopy(MassCompsNew, MassComps);
                        Fun = NewFun;
                    }
                    else
                    {
                        if (Math.Exp( - (NewFun - Fun) / Tk) > random())
                        {
                            ArrCopy(MassCompsNew, MassComps);
                            Fun = NewFun;
                        }
                    }
                }
            }
            // 程序计算完成
            return MassComps[0];
        }

• 测试代码

            double[] x = new double[1000];
            double[] y = new double[1000];
            double[] xfit = new double[201];
            double[] yfit = new double[201];
            double c = 5;
            double sigm = 3;
            for (int i = 0; i < 1000; i++)
            {
                x[i] = 0.01 * i;
                y[i] = Math.Exp(- (x[i] - c) * (x[i] - c) / (sigm * sigm));
                y[i] = y[i] + 0.1 * (random() - 0.5) * y[i];
            }
            for (int i = 0; i < 201; i++)
            {
                xfit[i] = x[400 + i];
                yfit[i] = y[400 + i];
            }
            double result = UnitPeakW.PeakPicker(xfit, yfit);
            MessageBox.Show(result.ToString());