Come adattare un gaussiano ai dati in matlab/ottava?

Question

Ho un set di dati di frequenza con picchi a cui ho bisogno di adattare una curva gaussiana e quindi ottenere la metà della larghezza massima da. La parte FWHM che posso fare, ho già un codice per questo, ma ho problemi a scrivere codice per adattarlo al gaussiano.

Qualcuno sa di funzioni che faranno questo per me o sarebbero in grado di indicarmi la giusta direzione? (Posso fare meno quadrati adatti per linee e polinomi ma non riesco a farlo funzionare per gaussiani)

Inoltre sarebbe utile se fosse compatibile con entrambi Octave e Matlab come ho Octave al momento ma don Non avrò accesso a Matlab fino alla prossima settimana.

Qualsiasi aiuto sarebbe molto apprezzato!

Answer 1

Montaggio direttamente una singola 1D gaussiana è un problema di montaggio non lineare. Troverete implementazioni pronti here, o here o here for 2D, o here (se avete la cassetta degli attrezzi statistiche) (avete sentito parlare di Google? :)

In ogni caso, ci potrebbe essere una soluzione più semplice. Se si sa per certo i dati y saranno ben descritti da una gaussiana, ed è ragionevolmente ben distribuite su tutto il x -range, è possibile linearizzare il problema (si tratta di equazioni, non dichiarazioni):

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

dove sono state fatte le sostituzioni

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

. Ora, risolvere per il sistema lineare Ax=b con (queste sono dichiarazioni 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;      

Il vettore x hai trovato in questo modo sarà pari

x == [P Q R] 

che è quindi necessario decodificare per trovare il significare μ e la deviazione standard σ:

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

Quale è possibile un controllo incrociato con x(3) == R (ci dovrebbe essere solo piccole differenze).

Answer 2

Forse questa è la cosa che stai cercando? Non sei sicuro di compatibilità: http://www.mathworks.com/matlabcentral/fileexchange/11733-gaussian-curve-fit

Dalla sua documentazione:

[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 3

ho avuto un problema simile. questo è stato il primo risultato su google, e alcuni degli script collegati qui hanno fatto il mio crash matlab.

finalmente ho trovato here che matlab ha una funzione integrata, che può adattarsi anche ai gaussiani.

sembrare che:

>> 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 4

ho trovato che il MATLAB "fit" funzione era lento, e utilizzati "lsqcurvefit" con una funzione gaussiana inline.Questo è per il montaggio di una FUNZIONE gaussiana, se vuoi semplicemente inserire i dati in una distribuzione normale, usa "normfit".

controllarlo

% % 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')