¿Cómo adaptar un gaussian a los datos en matlab/octava?

octave curve-fitting

Question

Tengo un conjunto de datos de frecuencia con picos a los que necesito ajustar una curva gaussiana y luego obtener la mitad máxima de ancho completo. La parte de FWHM que puedo hacer, ya tengo un código para eso, pero estoy teniendo problemas para escribir código que se ajuste al Gaussian.

¿Alguien sabe de alguna función que pueda hacer esto por mí o que sea capaz de orientarme en la dirección correcta? (Puedo hacer ajustes mínimos cuadrados para líneas y polinomios pero no puedo hacer que funcione para gaussianos)

También sería útil si fuera compatible con Octave y Matlab ya que tengo Octave en este momento pero no tengo acceso a Matlab hasta la semana próxima.

¡Cualquier ayuda sería muy apreciada!

Answer 1

tuve un problema similar este fue el primer resultado en google, y algunos de los scripts vinculados aquí hicieron que mi matlab fallara.

finalmente encontré aquí que el matlab tiene una función de ajuste incorporada, que también puede adaptarse a los gaussianos.

se ve así:

>> v=-30:30;
>> fit(v'', exp(-v.^2)'', ''gauss1'')

ans = 

   General model Gauss1:
   ans(x) =  a1*exp(-((x-b1)/c1)^2)
   Coefficients (with 95% confidence bounds):
      a1 =           1  (1, 1)
      b1 =  -8.489e-17  (-3.638e-12, 3.638e-12)
      c1 =           1  (1, 1)

Answer 2

Descubrí que la función de "ajuste" de MATLAB era lenta y usé "lsqcurvefit" con una función gaussiana en línea. Esto es para ajustar una FUNCIÓN Gaussiana, si solo quieres ajustar los datos a una distribución Normal, usa "normfit".

Revisalo

% % Generate synthetic data (for example) % % %

    nPoints = 200;  binSize = 1/nPoints ; 
    fauxMean = 47 ;fauxStd = 8;
    faux = fauxStd.*randn(1,nPoints) + fauxMean; % REPLACE WITH YOUR ACTUAL DATA
    xaxis = 1:length(faux) ;fauxData = histc(faux,xaxis);

    yourData = fauxData; % replace with your actual distribution
    xAxis = 1:length(yourData) ; 

    gausFun = @(hms,x) hms(1) .* exp (-(x-hms(2)).^2 ./ (2*hms(3)^2)) ; % Gaussian FUNCTION

% % Provide estimates for initial conditions (for lsqcurvefit) % % 

    height_est = max(fauxData)*rand ; mean_est = fauxMean*rand; std_est=fauxStd*rand;
    x0 = [height_est;mean_est; std_est]; % parameters need to be in a single variable

    options=optimset(''Display'',''off''); % avoid pesky messages from lsqcurvefit (optional)
    [params]=lsqcurvefit(gausFun,x0,xAxis,yourData,[],[],options); % meat and potatoes

    lsq_mean = params(2); lsq_std = params(3) ; % what you want

% % % Plot data with fit % % % 
    myFit = gausFun(params,xAxis);
    figure;hold on;plot(xAxis,yourData./sum(yourData),''k'');
    plot(xAxis,myFit./sum(myFit),''r'',''linewidth'',3) % normalization optional
    xlabel(''Value'');ylabel(''Probability'');legend(''Data'',''Fit'')

Answer 3

Quizás esto tiene lo que estás buscando? No estoy seguro acerca de la compatibilidad: http://www.mathworks.com/matlabcentral/fileexchange/11733-gaussian-curve-fit

De su documentación:

[sigma,mu,A]=mygaussfit(x,y) 
[sigma,mu,A]=mygaussfit(x,y,h)

this function is doing fit to the function 
y=A * exp( -(x-mu)^2 / (2*sigma^2) )

the fitting is been done by a polyfit 
the lan of the data.

h is the threshold which is the fraction 
from the maximum y height that the data 
is been taken from. 
h should be a number between 0-1. 
if h have not been taken it is set to be 0.2 
as default.

Answer 4

Montar un 1D Gaussian directamente es un problema de ajuste no lineal. Encontrarás implementaciones ya hechas aquí , o aquí , o aquí para 2D , o aquí (si tienes la caja de herramientas de estadísticas) (¿has oído hablar de Google? :)

De todos modos, podría haber una solución más simple. Si está seguro de que sus datos serán bien descritos por un gaussiano y están razonablemente bien distribuidos en todo su rango x , puede linealizar el problema (estas son ecuaciones, no declaraciones):

   y = 1/(σ·√(2π)) · exp( -½ ( (x-μ)/σ )² )
ln y = ln( 1/(σ·√(2π)) ) - ½ ( (x-μ)/σ )²
     = Px² + Qx + R         

donde las sustituciones

P = -1/(2σ²)
Q = +2μ/(2σ²)    
R = ln( 1/(σ·√(2π)) ) - ½(μ/σ)²

ha sido hecho. Ahora, resuelve para el sistema lineal Ax=b con (estas son declaraciones de Matlab):

% design matrix for least squares fit
xdata = xdata(:);
A = [xdata.^2,  xdata,  ones(size(xdata))]; 

% log of your data 
b = log(y(:));                  

% least-squares solution for x
x = A/b;                    

El vector x que encontraste de esta manera será igual

x == [P Q R]

que luego tiene que realizar una ingeniería inversa para encontrar la media μ y la desviación estándar σ:

mu    = -x(2)/x(1)/2;
sigma = sqrt( -1/2/x(1) );

Lo cual puedes verificar con x(3) == R (solo debería haber pequeñas diferencias).