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

[Stéphane Mottelet]

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
> users at lists.scilab.org
> 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20201109/c7db7a28/attachment.htm>