Regresión múltiple con funciones definidas por el usuario en Matlab

En ocasiones se nos presentan situaciones en las que se busca ajustar una función a un conjunto de datos, en los programas es habitual encontrarnos con funciones predefinidas para utilizarlas fácilmente, en Matlab se ofrecen algunas como:

‘poly1’    \rightarrow \; y=mx+c
‘power1’ \rightarrow \; y=ax^b
‘sin1’     \rightarrow \; y=a\cdot\sin(bx+c)

Estas configuraciones rápidas son muy útiles pero no nos sirven para todos los casos, en ocasiones se busca ajustar un conjunto de datos a otro tipo de función que se presupone que define el conjunto de datos y por tanto el ajuste será mejor. Teniendo de ejemplo una función como la de campo perjudicial planteada en la teoría revisada corregida, que se escribe del siguiente modo:

 

(1)   \begin{flalign*} L_{p,\text{perjudicial}}=y=10\log_{10}\left(\frac{Z}{p_0^2} \cdot\frac{4W}{A} \cdot e^{-\frac{13.82\left(\frac{r}{c}+t_0\right)\cdot \epsilon_L}{T}}\cdot C_L \right) \end{flalign*}

 

donde Z, p_0, W, A, c, t_0 y T son constantes, r es la variable independiente, \epsilon_L y C_L son los coeficientes de la función (los valores que se desean obtener) y L_{p,\text{perjudicial}} es el valor de salida de la función o y (en el caso del ejemplo, es la curva del campo perjudicial obtenido experimentalmente). Buscaremos, mediante regresión, que esta función devuelva unos valores similares a la curva ajustando los coeficientes definidos.

En Matlab la función que nos permite realizar una regresión mediante una función personalizada es fit, que para nuestro caso se escribirá como fit(x,y,fitType,Name,Value), como se puede observar se necesita crear un objeto de tipo fitType, en este objeto es donde estableceremos la ecuación para realizar la regresión, la variable dependiente (y), la variable independiente (r o x), las constantes (Z, p_0, W, A, c, t_0 y T) y los coeficientes (\epsilon_L y C_L) que en lenguaje Matlab es:

 

fP = fittype(...
       'Zimp/p0^2 * (4*W)/A * exp(-(13.82*(r/c+t0)*el/T)) * Cl',...   % Ecuación
       'dependent',{'y'},...                                          % Variable dependiente, salida de la ecuación
       'independent',{'r'},...                                        % Variable independiente, distancia        
       'problem', {'Zimp','p0','W','A','t0','T','c'},...              % Constantes
       'coefficients',{'el','Cl'});                                   % Coeficientes

 

Una vez definida nuestra regresión, para realizarla se deben introducir todas las variables y constantes y el objeto  fitType en la función fit, quedando:

 

[fitobject,gof,output] = fit(r,10.^(LPexp/10),fP,'problem',{Zimp,p0,W,A,t0,T,c});

 

La función fit puede devolver 3 variables:

  • fitobject es un objeto de tipo cfit similar a una estructura que contiene el valor de los coeficientes incluyendo un margen de confianza del 95%. También contiene la función/ecuación utilizada y el valor de los parámetros de entrada.
  • gof es un objeto de tipo struct que contiene los parámetros del ajuste de la regresión a los datos de entrada.
  • output es un objeto de tipo struct que contiene los detalles del proceso de la regresión, los residuos, etc.

La variable que nos interesa principalmente es la fitobject que nos dará los valores de los coeficientes a los que estábamos buscando valores.

Para comprender bien todo lo explicado (sé que es complejo de entender sin experimentar) incluyo el siguiente código que contiene todo lo necesario para realizar la regresión y obtener los resultados (al final de esta publicación incluyo la descarga de scripts para ejecutar directamente en Matlab, estos incluyen más información):

 

%% Ejemplo de regresión multivariable

% Variables
r       =   0:0.1:18;       % Vector de distancias (m)
LPexp   =   -0.23*r+75.28;  % Curva del campo perjudicial experimental (dB)
Zimp    =   413.48;         % Impedancia acústica del aire (rayls)
p0      =   2*10^-5;        % Presión acústica de referencia (pascales)
W       =   0.0026;         % Potencia acústica de la fuente (vatios)
A       =   66.3082;        % Absorción acústica equivalente (m^2)
t0      =   0.05;           % Tiempo de integración temporal (segundos)
T       =   1.1577;         % Tiempo de reverberación (segundos)
c       =   343.4;          % Velocidad del sonido en el aire (m/s)

%% Definición de la regresión (ecuación, variables, coeficientes, etc)
% La ecuación la introduciremos en lineal para facilitar el proceso de
% regresión, por lo que quitamos el logaritmo de la ecuación.
fP = fittype(...
    'Zimp/p0^2 * (4*W)/A * exp(-(13.82*(r/c+t0)*el/T)) * Cl',... % Ecuación
    'dependent',{'y'},...                                        % Variable dependiente, salida de la ecuación
    'independent',{'r'},...                                      % Variable independiente, distancia        
    'problem', {'Zimp','p0','W','A','t0','T','c'},...            % Constantes
    'coefficients',{'el','Cl'});                                 % Coeficientes

%% Regresión
% Tanto la variable dependiente como la independiente deben ser columnas:
r = r';
LPexp = LPexp';

% La variable de la curva del campo perjudicial experimental se convierte a
% lineal para que la comparación con la función introducida sea correcta,
% ya que antes hemos eliminado el logaritmo.
[fitobject,gof,output] = fit(r,10.^(LPexp/10),fP,'problem',{Zimp,p0,W,A,t0,T,c});

%% Variables obtenidas mediante la regresión
% Coeficientes
el = fitobject.el;
Cl = fitobject.Cl;
% Valor de ajuste R^2
R2 = gof.adjrsquare;

%% Curva obtenida mediante la regresión en decibelios
LPteo = 10*log10(Zimp/p0^2 * (4*W)/A * exp(-(13.82*(r/c+t0)*el/T)) * Cl);

%% Representación gráfica de ambas curvas
% Experimental
plot(r,LPexp,'r','LineWidth',2);
hold on;
% Teórico (regresión)
plot(r,LPteo,'b:','LineWidth',2);

 


Deja un comentario