Гауссовское соответствие в C#
В проекте, над которым я работаю, мне нужно получить гауссову подгонку из набора точек - нуждаются в среднем и дисперсии для некоторой обработки и, возможно, степени ошибки (или уровня точности), чтобы я мог выяснить, действительно ли набор точек имеют нормальное распределение.
Я нашел этот вопрос
но он ограничен только 3 баллами - тогда как мне нужно подгонка, которая может работать с любым количеством баллов.
Что мне нужно, это похоже на labview Gaussian Peak Fit
Я посмотрел на mathdotnet и aforge.net (используя оба в одном проекте), но ничего не нашел.
Кто-нибудь знает какие-либо C# или (легко конвертируемые) C/C++ или Java-решения?
В качестве альтернативы мне сказали, что следует использовать итеративный алгоритм - я мог бы реализовать его сам (если не слишком сложно по математике). Есть идеи о том, что я могу использовать? Я прочитал много статей (в Википедии и других, найденных через Google), но не нашел четкого указания на решение.
5 ответов
Я нашел хорошую реализацию в ImageJ, общедоступной программе обработки изображений, исходный код которой доступен здесь
Переоборудован в C# и адаптирован к моим потребностям.
Спасибо вам, ребята, за ваши ответы... не строго связанные с решением, которое я нашел, но каким-то образом я попал туда с вашей помощью:)
В Math.Net (NuGet), вы можете сделать:
var real_σ = 0.5;
var real_μ = 0;
//Define gaussian function
var gaussian = new Func<double, double, double, double>((σ, μ, x) =>
{
return Normal.PDF(μ, σ, x);
});
//Generate sample gaussian data
var data = Enumerable.Range(0, 41)
.Select(x => -2 + x * 0.1)
.Select(x => new { x, y = gaussian(real_σ, real_μ, x) });
var xs = data.Select(d => d.x).ToArray();
var ys = data.Select(d => d.y).ToArray();
var initialGuess_σ = 1;
var initialGuess_μ = 0;
var fit = Fit.Curve(xs, ys, gaussian, initialGuess_σ, initialGuess_μ);
//fit.Item1 = σ, fit.Item2 = μ
Здесь я показываю пример того, как вы можете подогнать произвольную функцию с произвольным количеством параметров с верхней / нижней границами для каждого отдельного параметра. Как и в случае с RexCardan, это делается с помощью библиотеки MathNet.
Это не очень быстро, но работает.
Чтобы соответствовать другой функции, измените
double gaussian(Vector<double> vectorArg)
метод. Если вы измените количество
vectorArg
s вам также необходимо настроить:
- Количество элементов в
lowerBound
,upperBound
иinitialGuess
вCurveFit
. - Измените количество параметров
return z => f(new DenseVector(new[] { parameters[0], parameters[1], parameters[2], parameters[3], parameters[4], parameters[5], z }));
- Измените количество параметров
t => f(new DenseVector(new[] { p[0], p[1], p[2], p[3], p[4], p[5], t }))
Пример кода для двойного гаусса
using MathNet.Numerics;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Linq;
static class GaussianFit
{
/// <summary>
/// Non-linear least square Gaussian curve fit to data.
/// </summary>
/// <param name="mu1">x position of the first Gaussian.</param>
/// <param name="mu2">x position of the second Gaussian.</param>
/// <returns>Array of the Gaussian profile.</returns>
public Func<double, double> CurveFit(double[] xData, double[] yData, double mu1, double mu2)
{
//Define gaussian function
double gaussian(Vector<double> vectorArg)
{
double x = vectorArg.Last();
double y =
vectorArg[0] * Normal.PDF(vectorArg[1], vectorArg[2], x)
+ vectorArg[3] * Normal.PDF(vectorArg[4], vectorArg[5], x);
return y;
}
var lowerBound = new DenseVector(new[] { 0, mu1 * 0.98, 0.05, 0, mu2 * 0.98, 0.05 });
var upperBound = new DenseVector(new[] { 1e10, mu1 * 1.02, 0.3, 1e10, mu2 * 1.02, 0.3 });
var initialGuess = new DenseVector(new[] { 1000, mu1, 0.2, 1000, mu2, 0.2 });
Func<double, double> fit = CurveFuncConstrained(
xData, yData, gaussian, lowerBound, upperBound, initialGuess
);
return fit;
}
/// <summary>
/// Non-linear least-squares fitting the points (x,y) to an arbitrary function y : x -> f(p0, p1, p2, x),
/// returning a function y' for the best fitting curve.
/// </summary>
public static Func<double, double> CurveFuncConstrained(
double[] x,
double[] y,
Func<Vector<double>, double> f,
Vector<double> lowerBound,
Vector<double> upperBound,
Vector<double> initialGuess,
double gradientTolerance = 1e-5,
double parameterTolerance = 1e-5,
double functionProgressTolerance = 1e-5,
int maxIterations = 1000
)
{
var parameters = CurveConstrained(
x, y, f,
lowerBound, upperBound, initialGuess,
gradientTolerance, parameterTolerance, functionProgressTolerance,
maxIterations
);
return z => f(new DenseVector(new[] { parameters[0], parameters[1], parameters[2], parameters[3], parameters[4], parameters[5], z }));
}
/// <summary>
/// Non-linear least-squares fitting the points (x,y) to an arbitrary function y : x -> f(p0, p1, p2, x),
/// returning its best fitting parameter p0, p1 and p2.
/// </summary>
public static Vector<double> CurveConstrained(
double[] x,
double[] y,
Func<Vector<double>, double> f,
Vector<double> lowerBound,
Vector<double> upperBound,
Vector<double> initialGuess,
double gradientTolerance = 1e-5,
double parameterTolerance = 1e-5,
double functionProgressTolerance = 1e-5,
int maxIterations = 1000
)
{
return FindMinimum.OfFunctionConstrained(
(p) => Distance.Euclidean(
Generate.Map(
x,
t => f(new DenseVector(new[] { p[0], p[1], p[2], p[3], p[4], p[5], t }))
),
y),
lowerBound,
upperBound,
initialGuess,
gradientTolerance,
parameterTolerance,
functionProgressTolerance,
maxIterations
);
}
пример
Чтобы уместить два гауссиана в положениях x 10 и 20:
Func<double, double> fit = GaussianFit.Curvefit(x_data, y_data, 10, 20);
Просто вычислите среднее значение и стандартное отклонение вашей выборки, это только два параметра гауссовского распределения.
Для "хорошего соответствия" вы можете сделать что-то вроде среднеквадратичной ошибки CDF.
Относительно ответа 1 Антонио 18 окт.
Я нашел хорошую реализацию в ImageJ , общедоступной программе обработки изображений, исходный код которой доступен здесь.
К сожалению, ссылка битая, однако на archive.org можно найти снимок:
https://imagej.nih.gov/ij/developer/source/ij/measure/CurveFitter.java.html
(Я бы опубликовал это как комментарий к ответу 1, но, как новый пользователь, мне, по-видимому, не разрешено это делать.)
Укко