[Scilab-users] Multiple regression on semi-log plot

classic Classic list List threaded Threaded
9 messages Options
arctica1963 arctica1963
Reply | Threaded
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

<http://mailinglists.scilab.org/file/t495709/Capture_Radial-averaged-power_spectra.jpg>



--
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
mottelet mottelet
Reply | Threaded
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)*a+(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

<https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/file/t495709/Capture_Radial-averaged-power_spectra.jpg> 



--
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
arctica1963 arctica1963
Reply | Threaded
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
<http://mailinglists.scilab.org/file/t495709/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
mottelet mottelet
Reply | Threaded
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)*a+(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
<https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/file/t495709/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].<a href="orghttps://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users" target="_blank">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
arctica1963 arctica1963
Reply | Threaded
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
mottelet mottelet
Reply | Threaded
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
Dang Ngoc Chan, Christophe Dang Ngoc Chan, Christophe
Reply | Threaded
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
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Re: Multiple regression on semi-log plot

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
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Re: Multiple regression on semi-log plot

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