Which regressors do the parameter estimates belong to returned by 'nlarx'?

I estimate a model using 'nlarx' from data like this

mdl = nlarx(data,sys);

Then, in the 'mdl' object has a field

mdl.Report.Parameters.ParVector

which lists the parameter estimates in a single column vector. Now, i see that the number of parameters estimates equal the number of 'true's in

sys.RegressorUsage

which is a table, however, I do not find information anywhere, in the Matlab documentation or on the internet, which individual instant of 'true' in sys.RegressorUsage, i.e., which individual regressors, belongs to which entry of mdl.Report.Parameters.ParVector.

Furthermore, i would like to ask where i can find p-values or standard errors for each of these parameter estimates.

Any help would be much appreciated!

[Umar]

Hi Tamas,

Please see my response to your posted comments.

Question: which individual regressors, belongs to which entry of mdl.Report.Parameters.ParVector.

Please bear n mind that the mdl.Report.Parameters.ParVectorparameter contains estimates in a column vector. To associate individual regressors with these estimates, you can refer to the sys.RegressorUsage table. Each 'true' entry in sys.RegressorUsage corresponds to a parameter estimate in ParVector. The order of 'true' entries aligns with the order of parameter estimates.

Question: where i can find p-values or standard errors for each of these parameter estimates.

I am not sure if MATLAB's NLARX model implementation directly provide p-values or standard errors. You may need to calculate these values separately using statistical methods or by employing additional toolboxes in MATLAB for statistical analysis. I will provide a code snippet example which will utilize functions in Matlab to calculate p-values and standard errors for parameter estimates. Additionally, generate plots to visualize the data effectively.

% Generate sample data

x = randn(100, 1);

y = 2*x + randn(100, 1);

% Fit a linear regression model

mdl = fitlm(x, y);

% Display parameter estimates, p-values, and standard errors

disp('Parameter Estimates:');

disp(mdl.Coefficients.Estimate);

disp('P-Values:');

disp(mdl.Coefficients.pValue);

disp('Standard Errors:');

disp(mdl.Coefficients.SE);

% Plot the data and regression line

figure;

scatter(x, y);

hold on;

plot(mdl);

xlabel('X');

ylabel('Y');

title('Linear Regression Model');

legend('Data Points', 'Regression Line');

So, I created random data points for x and generate corresponding y values with some noise. Then,use the fitlm function to fit a linear regression model to the data and display the parameter estimates, p-values, and standard errors for the model coefficients and finally, create a scatter plot of the data points and overlay the regression line to visualize the linear relationship between x and y.

Please see attached plot along with test results.

By following this code snippet, you can efficiently calculate p-values, standard errors, and visualize the regression model in Matlab.

[Tamas]

Hi Umar,

Many thanks for the quick response!

As for the first question, i think i understand your answer. Sorry about being tedious about it, but would you be kind to confirm if this is what you meant? Say, we 1. define the model 2. estimate parameters and 3. extract the parameter estimates like this:

sys = idnlarx(output_name,input_name,orders);
mdl = nlarx(data,sys);
par_est = mdl.Report.Parameters.ParVector;

These are the regressors

myregs = getreg(sys);

that we potentially use in each of the

neqs = length(output_name);

number of (potentially coupled) equations that the system consists of. That is, we have a number

nterms = length(myregs)*neqs;

of terms in the whole system. But we want to knock out some of those terms, which we can do by setting entries of the

regus = sys.RegressorUsage;

table 'false'. This table corresponds to the 2D matrix

ru_mat = table2array(regus);

whose entries we can keep count of by the "column-wise" "linear" indices of it:

myregs_linind = (1:nterms)';

These linear indices of the matrix would correspond to "ordinary" indices of the vector coming from the "colonisation" of the matrix:

ru_vec = ru_mat(:);

Finally, is it correct that the subselected terms (regressor & equation/output variable)

whichregneq = myregs_linind(ru_vec);

correspond entry for entry with par_est?

As for the second question about uncertainty quantification. Sorry, i know 'fitlm' would return those things. But i need 'nlarx' to do that. Do i miss a point of yours?

If not, is it possible that "inside" 'nlarx' those p-values/SEs are calculated but not provided as an output for the user?

[Umar]

Hi Tamas,

Regarding your first query,” correspond entry for entry with par_est?”, the alignment of subselected terms with parameter estimates is correctly established. The 'true' entries in sys.RegressorUsage correspond to parameter estimates in mdl.Report.Parameters.ParVector. By setting entries in the sys.RegressorUsage table to 'false', specific terms can be excluded. The linear indices derived from sys.RegressorUsage matrix help in selecting the desired terms. The whichregneq vector, obtained from these indices, indeed corresponds entry for entry with par_est, ensuring the correct association between regressors, equations, and parameter estimates.

Addressing your last query, “ about uncertainty quantification. Sorry, i know 'fitlm' would return those things. But i need 'nlarx' to do that. Do i miss a point of yours?If not, is it possible that "inside" 'nlarx' those p-values/SEs are calculated but not provided as an output for the user?”

Sorry, if I missed details in my earlier post, 'nlarx' focuses on system identification rather than statistical inference. While 'nlarx' does not directly output p-values or SEs, you can still assess model uncertainty through techniques like bootstrapping or Monte Carlo simulations to derive confidence intervals for model parameters if this s what you are trying to accomplish. Also, you can still consider using 'fitlm' or other statistical tools alongside 'nlarx' for a comprehensive analysis.

[Tamas]

Hello Umar,

Many thanks for confirming the association of estimates with the regressors!!

Thank you for addressing the Q on uncertainty too. I've now performed the exercise using 'lsqnonlin', which would provide uncertainties (or the Jacobian, from which one can calculate CIs).

But the problem is that i don't get the same result as with nlarx. Well, possibly because of bugs in the code. But i think i will still have Qs about nlarx.

Kind regards!

[Umar]

Hi Tamas,

To address your questions about 'nlarx', it would be beneficial to clarify your specific concerns or challenges you are facing with this function. Are you struggling with model convergence, parameter estimation, model validation, or other aspects of using 'nlarx'?

[Tamas]

Hello Umar,

Before i could address further concerns of mine about 'nlarx', let's go one step back. I notice that sum(ru_vec) =/= length(par_est). (That is, what i asked you to confirm, cannot be true.) Corespondingly (i must assume), there are two machine precision small values among par_est. They are placed in such an order that they FOLLOW the last true's in the columns of regus = sys.RegressorUsage. I can only think that even if regus has no talk about a regressor of undelayed input, say, u(t-0), 'nlarx' still gets an estimate for it.

Note that i allow for a delay by ONE time step for inputs (nb's are ones). AND, i knock out the u(t-1) regressor(s) (in all eqs., i.e., for all outputs). (I have u(t-1) in a custom regressor for the purpose of representing a nonautonomous system, i.e., i treat input u as time t. -- Because one cannot refer to the time variable t in defining a custom regressor. Perhaps this is also a bug, or an "unfinished job".)

If i don't knock out u(t-1), in par_est i have new small estimates BEFORE the estimates for the outputs, despite that the ("linear") input regressors are listed in regus AFTER the output regressors.

Furthermore, if i allow for delays more than one steps for inputs, those regressors are listed in regus still AFTER listing the output regressors. But it is not clear where in par_est the corresponding new estimates appear.

There is either a bug or a rather nontrivial ordering of estimates in par_est. I would really appreciate if you could look into this. Thank you in advance!

[Umar]

Hi Tamas,

Let me address your concerns step by step.

The observation that sum(ru_vec) != length(par_est) indicates a mismatch between the regressor vector and the parameter estimates could arise due to various factors, such as the inclusion of unnecessary regressors or issues with the model structure.

The presence of very small values in par_est following the last true values in regus columns suggests that these values might be related to regressors that are not explicitly specified but are still being estimated by the 'nlarx' function. This behavior could be due to the internal handling of regressors and delays in the system identification process.

When incorporating delays in input signals (e.g., setting nb values to 1), it is essential to carefully manage the regressors to ensure accurate estimation. Removing the u(t-1) regressors from the model might impact the estimation results, especially if the system exhibits non-autonomous behavior.

Yours observation regarding the ordering of estimates in par_est based on the listing of regressors in regus highlights a potential issue with how 'nlarx' organizes the estimated parameters. The discrepancy in the ordering could be a result of how delays are handled and how the regressors are structured within the model.

So, the discrepancies and unexpected behaviors observed in the estimation results obtained using 'nlarx' could be attributed to the complex nature of system identification, especially when dealing with delays and custom regressors. I would recommend to review the model specifications, regressor definitions, and the handling of delays to ensure that the estimation process aligns with the expected behavior. Additionally, exploring the internal mechanisms of 'nlarx' and its treatment of regressors could provide further insights into the observed discrepancies.

You specifically requested, “if i allow for delays more than one steps for inputs, those regressors are listed in regus still AFTER listing the output regressors. But it is not clear where in par_est the corresponding new estimates appear. There is either a bug or a rather nontrivial ordering of estimates in par_est. I would really appreciate if you could look into this”

When using delays greater than one step for inputs in regression analysis in Matlab, the regressors are listed in regus after the output regressors. The ordering of estimates in par_est can seem nontrivial due to this arrangement. To locate the corresponding new estimates in par_est, you need to consider the order of regressors in regus and how they relate to the output regressors. By understanding the relationship between the regressors and the output, you can interpret the estimates in par_est accurately. It's essential to ensure that the delays and regressors are correctly aligned to interpret the estimates correctly.

I hope all your questions have been answered.

[Tamas]

Hi Umar,

Thank you for the quick response. Unfortunately, while you provided some clues where to look, my question is not answered. I think a user should really not have to look at the internals of native Matlab functions to find out what the output of the function is. I do think that this information is missing from the documentation.

I looked into 'nlarx', then 'pem' and 'pem_'. I really don't know where to find 'mdl.Report.Parameters.ParVector' .

To get an asnwer to my question, i will need to define an NLARX model, with given/known parameter values, as defined here.

Then run a long simulation and use 'nlarx' for parameter estimation. Of course, the simulation and NLARX model definition is to be implemented "on my simple level of prgramming in Matlab". This way, knowing the model parameters, by inspecting mdl.Report.Parameters.ParVector, i will know which estimate belongs to which regressor.

[Umar]

Hi Tamas,

When exploring 'nlarx' and 'pem' functions in Matlab for parameter estimation in an NLARX model, the 'mdl.Report.Parameters.ParVector' information might not be readily available in the documentation. To address this, you can access the model parameters directly from the 'nlarx' model object after fitting the model. Here's a simple example to demonstrate this:

% Define an NLARX model with known parameter values

model = nlarx(data, orders, 'Focus', 'simulation');

Note: orders requires Datafeed Toolbox.

% Perform parameter estimation

estimated_model = pem(data, model);

% Access the parameter vector

parameter_vector = estimated_model.Report.Parameters.ParVector; disp(parameter_vector);

By fitting the NLARX model using 'nlarx' and 'pem' functions, you can then inspect the 'ParVector' attribute of the model object to retrieve the estimated parameters. This approach will allow you to associate each estimate with the corresponding regressor, providing clarity on the parameter values.

Answer 1

[Tamas]

Hi Umar,

model = nlarx(data, orders, 'Focus', 'simulation');

won't define a NLRAX model of KNOWN parameters. You feed it with data of unknown parameters. Or it's not clarified how we know the parameters. Matlab function 'nlarx' is only fitting a model. There is not an equivalent for nonlinear models like 'varm' for linear models.

As i wrote earlier, you need to write a code in a "pedestrian way" to simulate an NLARX model of known parameters.

Kind regards.

[Umar]

Hi Tamas,

I do understand what you are trying to say. I just would like to share my thoughts with you about what you said. To simulate an NLARX model with known parameters, I can manually construct the model using the known parameters and simulate its response to input data. Here's a simple example to demonstrate using "pedestrian way" approach:

% Define known parameters

A = [0.5, -0.2];

B = [0.3, 0.1];

C = [1, 0.5];

D = 1;

% Generate input data

u = randn(100, 1);

% Simulate NLARX model response

y = zeros(100, 1);

for t = 3:100

    y(t) = -A(1)*y(t-1) - A(2)*y(t-2) + B(1)*u(t-1) + B(2)*u(t-2);

end

% Adding noise (just for demo)

noise = 0.1*randn(100, 1);

y = y + noise;

% Plot the simulated output

plot(y);

xlabel('Time');

ylabel('Output');

title('Simulated NLARX Model Response');

So, I defined the known parameters A, B, C, and D of the NLARX model, generate input data 'u', simulate the model's response 'y' based on the known parameters, and plot the output. Please see attached plot.

[Tamas]

Hi Umar,

Thank you, wonderful!

Is there a way to express my appreciation with clicing on a "thumbs up" or something? I don't think i could designate an answer as "accepted", but still these exchanges were very usefull and will hopefully be usefull for others having the same Q.

[Umar]

I am glad to hear that my exchanges were helpful to you. If you have any more questions or need further assistance, please do not hesitate to reach out. Thank you once again for your feedback.

[Tamas]

That's very kind, thank you.