# [Scilab-users] Multiple regression on semi-log plot Classic List Threaded 9 messages Open this post in threaded view
|

## [Scilab-users] Multiple regression on semi-log plot Hello, I am looking to determine multiple regression lines from a single power spectral dataset (log power vs radial wavenumber), and was wondering if it is feasible in Scilab to compute something similar to the attached plot? I did locate a Matlab code for finding a turning point in a plot and perhaps this may be a route to find the change in curvature? clear clf() //to find the "knee" dt = 0.01; t = 0:dt:1; y = exp(-10*t); //Compute first and second derivatives by finite differencing (centred) //yp = nan(size(y)); //ypp = nan(size(y)); yp=%nan*y; // similar behaviour to Matlab nan(size(y)) ypp=%nan*y; yp(2:\$-1) = (y(3:\$) - y(1:\$-2))/(2*dt); // first derivative ypp(2:\$-1) = (y(3:\$) + y(1:\$-2) - 2*y(2:\$-1)) / (dt^2); // second derivative //Compute the curvature k = abs(ypp) ./ (1 + yp.^2).^1.5 //Find the maximum curvature point and plot it on the curve [kmax, idx] = max(k); plot(t, y, 'b', t(idx), y(idx), 'ro') //plot the curvature figure() plot(t, k) Any pointers or ideas would be welcome. Thanks Lester -- Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Multiple regression on semi-log plot Hi, This is an easy task that can be done by fitting a piecewise-affine function like this: y = a*(x-theta)+phi, for x < theta y = b*(x-theta)+phi, for x >= theta Here is an example : ```function y=fun(x, param) a = param(1); b = param(2); theta = param(3); phi = param(4); c = (x=theta)*b; y = c.*(x-theta)+phi; endfunction function r=residual(param, x, y, f) r = sum((y-f(x,param)).^2); endfunction function [r, dr, ind]=costf(p, ind, x, y, f) r = residual(p,x,y,f); dr = numderivative(list(residual,x,y,fun),p); endfunction x=linspace(0,4,50); theta = 1.5; phi = 1; a = grand(1,1,"unf",-2,2) b = grand(1,1,"unf",-2,2) y = fun(x,[a,b,theta,phi]) + rand(x,"normal")/5; [ropt,popt] = optim(list(costf,x,y,fun),[0,0,mean(x),mean(y)]); clf plot(x,y,'o',x,fun(x,popt),popt(3),popt(4),'xr')``` S. Le 09/11/2020 à 17:07, arctica1963 a écrit : ```Hello, I am looking to determine multiple regression lines from a single power spectral dataset (log power vs radial wavenumber), and was wondering if it is feasible in Scilab to compute something similar to the attached plot? I did locate a Matlab code for finding a turning point in a plot and perhaps this may be a route to find the change in curvature? clear clf() //to find the "knee" dt = 0.01; t = 0:dt:1; y = exp(-10*t); //Compute first and second derivatives by finite differencing (centred) //yp = nan(size(y)); //ypp = nan(size(y)); yp=%nan*y; // similar behaviour to Matlab nan(size(y)) ypp=%nan*y; yp(2:\$-1) = (y(3:\$) - y(1:\$-2))/(2*dt); // first derivative ypp(2:\$-1) = (y(3:\$) + y(1:\$-2) - 2*y(2:\$-1)) / (dt^2); // second derivative //Compute the curvature k = abs(ypp) ./ (1 + yp.^2).^1.5 //Find the maximum curvature point and plot it on the curve [kmax, idx] = max(k); plot(t, y, 'b', t(idx), y(idx), 'ro') //plot the curvature figure() plot(t, k) Any pointers or ideas would be welcome. Thanks Lester -- Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html _______________________________________________ users mailing list [hidden email] https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users ``` ```-- Stéphane Mottelet Ingénieur de recherche EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet ``` _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Multiple regression on semi-log plot Hello, Thanks for the idea and suggestions. Not too sure how to apply it, if you could give some pointers on the attached data and code. The ultimate idea is to get the slopes of the straight line segments. Many thanks, Lester clear clf() // Read data - wavelength (in km)), power, 1 standard deviation // Unknown data length; 3 columns -default space delimited // PSD_wavelength.dat from GMT grdfft radially averaged power spectra data = read('PSD_wavelength.dat',-1,3); wavelength = data(:,1); power = data(:,2); std_dev1 = data(:,3); ln_power = log(power); wavenumber = 1./wavelength; f=gcf(); //plot(wavenumber, ln_power) scatter(wavenumber, ln_power,'marker','.') a=gce().children; a.mark_mode = "on" a.mark_style = 0 a.mark_size_unit = "point" a.mark_size=3 xlabel ('wavenumber k (km-1)') ylabel ('Log (Power)') PSD_wavelength.dat   -- Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Multiple regression on semi-log plot Hello, You just have to replace "x" by "wavelength" and "y" by "ln_power". The slopes are the first two component of the optimal vector "popt" : ```clear clf() // Read data - wavelength (in km)), power, 1 standard deviation // Unknown data length; 3 columns -default space delimited // PSD_wavelength.dat from GMT grdfft radially averaged power spectra data = read('PSD_wavelength.dat',-1,3); wavelength = data(:,1); power = data(:,2); std_dev1 = data(:,3); ln_power = log(power); wavenumber = 1./wavelength; f=gcf(); //plot(wavenumber, ln_power) scatter(wavenumber, ln_power,'marker','.') a=gce().children; a.mark_mode = "on" a.mark_style = 0 a.mark_size_unit = "point" a.mark_size=3 xlabel ('wavenumber k (km-1)') ylabel ('Log (Power)') function y=fun(x, param) a = param(1); b = param(2); theta = param(3); phi = param(4); c = (x=theta)*b; y = c.*(x-theta)+phi; endfunction function r=residual(param, x, y, f) r = sum((y-f(x,param)).^2); endfunction function [r, dr, ind]=costf(p, ind, x, y, f) r = residual(p,x,y,f); dr = numderivative(list(residual,x,y,fun),p); endfunction [ropt,popt] = optim(list(costf,wavenumber,ln_power,fun),[0,0,mean(wavenumber),mean(ln_power)]); plot(wavenumber,fun(wavenumber,popt),'r',popt(3),popt(4),'xr')``` arctica1963 <[hidden email]> a écrit : Hello, Thanks for the idea and suggestions. Not too sure how to apply it, if you could give some pointers on the attached data and code. The ultimate idea is to get the slopes of the straight line segments. Many thanks, Lester clear clf() // Read data - wavelength (in km)), power, 1 standard deviation // Unknown data length; 3 columns -default space delimited // PSD_wavelength.dat from GMT grdfft radially averaged power spectra data = read('PSD_wavelength.dat',-1,3); wavelength = data(:,1); power = data(:,2); std_dev1 = data(:,3); ln_power = log(power); wavenumber = 1./wavelength; f=gcf(); //plot(wavenumber, ln_power) scatter(wavenumber, ln_power,'marker','.') a=gce().children; a.mark_mode = "on" a.mark_style = 0 a.mark_size_unit = "point" a.mark_size=3 xlabel ('wavenumber k (km-1)') ylabel ('Log (Power)') PSD_wavelength.dat -- Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html _______________________________________________ users mailing list [hidden email].orghttps://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Multiple regression on semi-log plot Hello, Thanks for showing how it works, always good to get a proper understanding. Just wondering if the code can be adapted to find multiple straight line segments? Looking at my original test data, it appears to have "kinks" at wavenumber (0.05, 0.12 and 0.16), the last one is where the x plots. I am guessing that it needs some threshold or tolerance control as the angles at the kink points may be small or larger. Many thanks for the guidance. Lester -- Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: Multiple regression on semi-log plot Hi, It is possible to fit a more general piecewise-affine model, but you will have to know in advance the number of "kinks" S. Le 16/11/2020 à 09:49, arctica1963 a écrit : > Hello, > > Thanks for showing how it works, always good to get a proper understanding. > > Just wondering if the code can be adapted to find multiple straight line > segments? Looking at my original test data, it appears to have "kinks" at > wavenumber (0.05, 0.12 and 0.16), the last one is where the x plots. I am > guessing that it needs some threshold or tolerance control as the angles at > the kink points may be small or larger. > > Many thanks for the guidance. > Lester > > > > -- > Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html> _______________________________________________ > users mailing list > [hidden email] > https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users-- Stéphane Mottelet Ingénieur de recherche EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users
Open this post in threaded view
|

## Re: {EXT} Re: Multiple regression on semi-log plot Hello, > De : users <[hidden email]> De la part de Stéphane Mottelet > Envoyé : lundi 16 novembre 2020 09:56 > > It is possible to fit a more general piecewise-affine model, > but you will have to know in advance the number of "kinks" Well, I think it is possible to automatically adjust the number of segments: 1. Perform the regression with only n = 1 line segment. 2. Increase the number of segments by 1, n = n + 1. Perform a new regression. You can detect the furthest point to know where to separate the segments. 3. Compare the quadratic errors between the current and the preceding regression. If the difference is small enough then end else go back to 2. Regards -- Christophe Dang Ngoc Chan Mechanical calculation engineer General This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden. _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users Hi, Le 16/11/2020 à 10:26, Dang Ngoc Chan, Christophe a écrit : > Hello, > >> De : users <[hidden email]> De la part de Stéphane Mottelet >> Envoyé : lundi 16 novembre 2020 09:56 >> >> It is possible to fit a more general piecewise-affine model, >> but you will have to know in advance the number of "kinks" > Well, I think it is possible to automatically adjust the number of segments: > > 1. Perform the regression with only n = 1 line segment. > > 2. Increase the number of segments by 1, n = n + 1. > Perform a new regression. > You can detect the furthest point to know where to separate the segments. > > 3. Compare the quadratic errors between the current and the preceding regression. > If the difference is small enough then end That's the more difficult part. Considering the provided data, which will never fit well this kind of model, the quadratic error will decrease quite regularly. If the real concern is about finding the points with maximum curvature, a non parametric model should be more adequate (Local regression, Kernel regression, splines, or wathever). S. > else go back to 2. > > Regards > > -- > Christophe Dang Ngoc Chan > Mechanical calculation engineer > > General > This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden. > _______________________________________________ > users mailing list > [hidden email] > https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users-- Stéphane Mottelet Ingénieur de recherche EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users Hello, Just a quick thought. But could a moving window over the wavenumber provide a simpler route to locating the kink points for computing the slope segments? Perhaps a window size of say 0.05 and shifted 0.025 with the existing methodology you outlined Stephane. Lester -- Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users