Original file (SVG file, nominally 478 × 364 pixels, file size: 23 KB)
This is a file from the Wikimedia Commons. Information from its description page there is shown below. Commons is a freely licensed media file repository. You can help. |
Summary
DescriptionRegression pic gaussien dissymetrique bruite.svg |
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
This media was created with Scilab, a free open-source software. Here is a listing of the Scilab source used to create this file. |
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
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.http://www.gnu.org/copyleft/fdl.htmlGFDLGNU Free Documentation Licensetruetrue |
- 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.
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/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 10:50, 20 November 2012 | 478 × 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... |
File usage
The following page uses this file:
Global file usage
The following other wikis use this file:
- Usage on fr.wikipedia.org
Metadata
This file contains additional information, probably added from the digital camera or scanner used to create or digitize it.
If the file has been modified from its original state, some details may not fully reflect the modified file.
Short title | Regression sur un pic gaussien dissymétrique bruité |
---|---|
Image title | Creator: GL2PS 1.3.2, (C) 1999-2006 Christophe Geuzaine (geuz@geuz.org)
For: Scilab CreationDate: Tue Nov 20 11:41:31 2012 |
Width | 478.01855 |
Height | 364.31445 |