Jump to content

File:Regression pic gaussien dissymetrique bruite.svg

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikipedia, the free encyclopedia

Original file (SVG file, nominally 478 × 364 pixels, file size: 23 KB)

Summary

Description
English: Profile fitting of an asymmetrical noisy gaussian peak, by Gauss-Newton regression. Created with Scilab, processed with Inkscape.
Français : Ajustement de profil d'un pic gaussien dissymétrique bruité, par régression en utilisant l'algorithme de Gauss-Newton. Créé par Scilab, modifié par Inkscape.
Date
Source Own work
Author Cdang
Other versions Animated GIF: File:Regression pic assymetrique.gif. Savitzky-Golay smoothing and derivation of raw data: File:Savitzky-golay pic gaussien dissymetrique bruite.svg

Scilab source

English: English version by default.
Français : Version française, si les préférences de votre compte sont réglées (voir Special:Preferences).

File common_functions.sce.

function [y] = gauss_dissym(A, x)
    // generates a dissymetrical gaussian peak
    // inputs: A(1): peak position
    //     A(2): height of the peak
    //     A(3): width on the left side
    //     A(4): width on the right side
    //     x: vector of numbers
    // oututs: y: vector of numbers
    index = (x < A(1)); // vector of %t at the left, %f at the right
    y = zeros(x); // initialisation
    y(index) = A(2)*exp(-(x(index) - A(1)).^2/A(3)); // left profile
    y(~index) = A(2)*exp(-(x(~index) - A(1)).^2/A(4)); // right profile
endfunction

function [B] = approximative_linearisation(f, A, x, deltaA)
    // computes an approximation of the partial derivatives of f
    // with respect to its parameters A(i)
    // f: function of x which parameters are described by the vector A
    // as following: f = (A,x)
    // A: vector of parameters (numbers)
    // x: vector of numbers
    sizex = size(x, '*');
    sizeA = size(A, '*');
    B = zeros(sizex, sizeA);
    for i = 1:sizeA // only the parameter A(i) is modified
        Aleft = [A(1:(i-1)), A(i) - deltaA, A(i+1:$)];
        Aright = [A(1:(i-1)), A(i) + deltaA, A(i+1:$)];
        B(:,i) = (f(Aright, x) - f(Aleft, x))/2/deltaA; // df/dA(i)
    end
endfunction

File asymmetric_gauss_peak_generator.sce: generates the file noisy_asym_gauss_peak.txt.

// **********
// Constants and initialisation
// **********

clear;
clf;
chdir('mypath/');

// parameters of the noisy curve
paramgauss(1) = 0; // peak position
paramgauss(2) = 10; // height of the peak
paramgauss(3) = 1; // width on the left side of the peak
paramgauss(4) = 0.3; // width on the right side of the peak
var = 0.5; // variance of the normal law (noise)
nbpts = 200 // number of points
halfwidth = 3*max(paramgauss(3), paramgauss(4)); // to determine the x range
step = 2*halfwidth/nbpts;

// **********
// Functions
// **********

exec('common_functions.sce', -1)

// **********
// Main program
// **********

// Data generation

for i = 1:nbpts
    x = step*i - halfwidth;
    Xinit(i) = x;
    noise(i) = var*rand(1, 'normal');
end

Ybasis = gauss_dissym(paramgauss, Xinit);
Yinit = Ybasis + noise;

// Data saving

//plot(Xinit, Ybasis, "r")
//plot(Xinit, Yinit, "b")

write ('noisy_asym_gauss_peak.txt', [Xinit, Yinit])

File asymmetric_peak_regression.sce: processes the file noisy_asym_gauss_peak.txt and determines the parameters of the model function by regression. The initial parameters are retrieved from the Savitzky-Golay smooting and derivation, see File:Savitzky-golay pic gaussien dissymetrique bruite.svg.

// **********
// Constants and initialisation
// **********

clear;
clf;

chdir('mypath/')

// Newton-Raphson parameters
precision = 1e-7; // stop condition
itermax = 30; // idem

// Precision of approximated derivation
epsilon = 1e-6;

// **********
// Functions
// **********

exec('common_functions.sce', -1)

function [e] = res(Yexp, Ycal)
    e = sqrt(sum((Yexp-Ycal).^2));
endfunction

function [A, R] = gaussnewton(f, X, Yexp, A0, imax, epsilon)
    // A: set of parameters of the model optimised by regression (vector)
    // R: list of the quality factors for the regression (vector)
    // X: explanatory variable (vector)
    // Yexp : response variable, measured values (vector)
    // A0 : paramètres d'initialisation du modèle (vector)
    // epsilon : stop value (scalar)
    k = 1; // initial damping factor, <=1,
    // avoids divergence 
    n = size(X, '*');
    e0 = sqrt(sum(Yexp.^2)); // normalisation of the quality factor
    Ycal = f(A0, X); // nitial model
    R(1) = res(Yexp, Ycal)/e0; // initial quality factor
    disp('i = 1 ; R = '+string(R(1))) // display of initial param
    i = 1;
    B = A0;
    flag = %t
    while (i < imax) & flag // tests the global convergence
        i = i+1;
        deltay = Yexp - Ycal;
        J = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        tJ = J'; // transposed
        deltap0 = inv((tJ*J))*tJ*deltay;
        flag2 = %t // for the 1st execution
        while flag2 & (k>0.1)
            deltap = k*deltap0;
            Bnew = B + deltap';
            Ycal = f(Bnew, X);
            R(i) = res(Yexp, Ycal)/e0;
            flag2 = (R(i) >= R(i-1)) // true if it diverges
            if flag2 then k = k*0.75; // increase the damping on divergence
            else k0 = k; // to display the value
                k = (1 + k)/2; // reduces the damping on convergence
            end
        end
        B = Bnew;
        flag = abs(R(i-1) - R(i)) > epsilon;
        disp('i = '+string(i)+' ; R = '+string(R(i)))
    end
    A = B;
endfunction

// **********
// Main program
// **********

// Data reading
data = read('noisy_asym_gauss_peak.txt',-1,2);

// Data characteristics
Xdef = data(:,1);
Ydef = data(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];

// Regression
tic();
[Aopt, Rnr] =...
    gaussnewton(gauss_dissym, Xdef, Ydef,...
    Ainit, itermax, precision)
t = toc();

// Calculated curve

Yopt = gauss_dissym(Aopt, Xdef);

// Display

print(%io(2),Ainit)
print(%io(2),Aopt)
print(%io(2),t)

clf

subplot(2,1,1)
plot(Xdef, Ydef, "-b")
plot(Xdef, Yopt, "-r")

subplot(2,1,2)
plot(Rnr)

It is possible to compute the pseudo-inverse matrix using the Scilab function pinv. The pseudo-inversion is performed by singular values decomposition. The performance may vary according to the data structure. The lines

        J = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        tJ = J'; // transposed
        deltap0 = inv((tJ*J))*tJ*deltay;

must be replaced by

        J = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        deltap0 = pinv(J)*deltay;

It is also possible to use the Scilab function leastsq: the code is more compact, but becomes more "black box". Additionally, it is not possible to display the progression of the fiability factor R, and it is slower (141 ms vs 16 ms, i.e. 8.8 times slower).

// **********
// Constants and initialisation
// **********

clear;
clf;

chdir('mypath/')

// **********
// Functions
// **********

exec('common_functions.sce', -1)

function [e] = res(A, X, Yexp)
    Ycal = gauss_dissym(A, X);
    e = sqrt(sum((Yexp-Ycal).^2));
endfunction

// **********
// Main program
// **********

// Data reading
data = read('noisy_asym_gauss_peak.txt',-1,2);

// Data characteristics
Xdef = data(:,1);
Ydef = data(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];

// Regression
tic();
[Rnr, Aopt] =...
    leastsq(list(res, Xdef, Ydef), Ainit)
t = toc();

// Calculated curve

Yopt = gauss_dissym(Aopt, Xdef);

// Display

print(%io(2),Ainit)
print(%io(2),Aopt)
print(%io(2),t)

plot(Xdef, Ydef, "-b")
plot(Xdef, Yopt, "-r")

Licensing

I, the copyright holder of this work, hereby publish it under the following licenses:
GNU head Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled GNU Free Documentation License.
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported, 2.5 Generic, 2.0 Generic and 1.0 Generic license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
You may select the license of your choice.

Captions

Add a one-line explanation of what this file represents

Items portrayed in this file

depicts

20 November 2012

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current10:50, 20 November 2012Thumbnail for version as of 10:50, 20 November 2012478 × 364 (23 KB)Cdang{{Information |Description ={{en|1=Profile fitting of an asymmetrical noisy gaussian peak, by Newton-Raphson regression. Created with Scilab, processed with Inkscape.}} {{fr|1=Ajustement de profil d'un pis gaussien dissymétrique bruité, par régre...

The following page uses this file:

Global file usage

The following other wikis use this file:

Metadata