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

[stephane.mottelet at utc.fr]

  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" :
clearclf()// 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 = 0a.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 <arctica1963 at gmail.com> 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
> users at lists.scilab.orghttps://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20201115/7fc61a32/attachment.htm>