Integral2 seems to substitute non-scalar values of variable into integrand. Why?

Dear community,
I am struggling with an error that my code produces. I am trying to do a maximum-likelihood estimation (MLE). I am importing CSV data; each row is an observation and an observation is four-dimensional with column names "bad_bought", "good_bought", "bad_notbought", "good_notbought".

Each observation is assumed to have latent variables and that are drawn from a joint-normal distribution with mean and variance and and correlation . An observation’s numeric values for the four columns are given by
"bad_bought":
"good_bought":
"bad_notbought":
"good_notbought":
where the error terms are iid, normal, and have variance .

I understand that I can calculate the joint distribution of the four answers (for given ) directly. However, I want to use numerical integration, as in the code below, because the model will eventually become more complicated necessitating it.

I receive two type of error messages:
First, and this is copy-pasted below, MATLAB tells me my mvnpdf arguments do not have the right dimensions. I do not understand that. [t1, t2], mu_vec are both of dimensions 1×2, i.e., they have the same number of columns. This error disappears if I replace the first argument of mvnpdf with [20, 20] (or some other numbers). I do not understand why; should integral2 not substitute scalar values for t1 and t2?

Second, and this appears if I replace [t1,t2] with [20,20] in the argument of mvnpdf, I receive the error message telling me I cannot use "*" to multiply the normpdf values together. This suggests to me that again t1 and t2 are not scalars.

Could you help? Thank you!

% Load necessary data
data = readtable(‘df_restricted.csv’);

% Select relevant columns and convert to matrix
data_matrix = table2array(data(:, {‘bad_bought’, ‘good_bought’, ‘bad_notbought’, ‘good_notbought’}));

% Define the negative log-likelihood function
function neg_log_lik = neg_log_lik_real(params, data_matrix)

% Extract parameters
mu1 = params(1);
mu2 = params(2);
v1 = params(3);
v2 = params(4);
rho = params(5);
v = params(6);

% Mean vector
mu_vec = [mu1, mu2];
disp(size(mu_vec))

% Covariance matrix
cov_matrix = [v1, rho * sqrt(v1 * v2); rho * sqrt(v1 * v2), v2];
disp(size(cov_matrix))
log_sum = 0;

% Loop through each data point
for i = 1:size(data_matrix, 1)

% Define the integrand for 2-dimensional integration
f = @(t1, t2) ( …
mvnpdf([t1, t2], mu_vec, cov_matrix) …
* normpdf(data_matrix(i, 1) – (t1 – t2), 0, sqrt(v)) …
* normpdf(data_matrix(i, 2) – (t1 + t2), 0, sqrt(v)) …
* normpdf(data_matrix(i, 3) – (t1 + t2), 0, sqrt(v)) …
* normpdf(data_matrix(i, 4) – (t1 – t2), 0, sqrt(v)) …
);
% Before integral2 call

% Perform 2-dimensional numerical integration using integral2
pdf_value = integral2(f, -Inf, Inf, -Inf, Inf, ‘AbsTol’, 1e-6, ‘RelTol’, 1e-6);

% Take log of the result and accumulate
log_value = log(pdf_value);
log_sum = log_sum + log_value;
end

% Return negative log-likelihood
neg_log_lik = -log_sum;
end

% Set initial parameters
initial_params = [56, 11, 119, 131, -0.36, 419];

% Define the function handle for optimization
neg_log_lik_handle = @(params) neg_log_lik_real(params, data_matrix);

% Use fminunc to minimize the negative log-likelihood
options = optimoptions(‘fminunc’, ‘Algorithm’, ‘quasi-newton’, ‘MaxIterations’, 500, ‘TolFun’, 1e-8);
[estimated_params, fval, exitflag, output] = fminunc(neg_log_lik_handle, initial_params, options);

% Display the results
disp(‘Estimated parameters:’);
disp(estimated_params);

disp(‘Final log-likelihood value:’);
disp(-fval);

First error:

Error using mvnpdf (line 67)
X and MU must have the same number of columns.

Error in
Structural_Estimation>@(t1,t2)(mvnpdf([t1,t2],mu_vec,cov_matrix)*normpdf(data_matrix(i,1)-(t1-t2),0,sqrt(v))*normpdf(data_matrix(i,2)-(t1+t2),0,sqrt(v))*normpdf(data_matrix(i,3)-(t1+t2),0,sqrt(v))*normpdf(data_matrix(i,4)-(t1-t2),0,sqrt(v)))
(line 34)
mvnpdf([t1, t2], mu_vec, cov_matrix) …

First error:

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), …

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, …

Error in integralCalc (line 113)
[q,errbnd] = vadapt(@minusInfToInfInvTransform,interval, …

Error in
integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions)
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( …

Error in
integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( …

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, …

Error in integralCalc (line 113)
[q,errbnd] = vadapt(@minusInfToInfInvTransform,interval, …

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in Structural_Estimation>neg_log_lik_real (line 43)
pdf_value = integral2(f, -Inf, Inf, -Inf, Inf, ‘AbsTol’, 1e-6, ‘RelTol’, 1e-6);

Error in Structural_Estimation>@(params)neg_log_lik_real(params,data_matrix) (line 58)
neg_log_lik_handle = @(params) neg_log_lik_real(params, data_matrix);

Error in fminunc (line 233)
f = feval(funfcn{3},x,varargin{:});

Error in Structural_Estimation (line 62)
[estimated_params, fval, exitflag, output] = fminunc(neg_log_lik_handle, initial_params, options);

Caused by:
Failure in initial objective function evaluation. FMINUNC cannot continue.
—-
Second error:
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix
matches the number of rows in the second matrix. To operate on each element of the matrix
individually, use TIMES (.*) for elementwise multiplication.

Error in Structural_Estimation>@(t1,t2)(mvnpdf([20,20],mu_vec,cov_matrix)*normpdf(data_matrix(i,1)-(t1-t2),0,sqrt(v))*normpdf(data_matrix(i,2)-(t1+t2),0,sqrt(v))*normpdf(data_matrix(i,3)-(t1+t2),0,sqrt(v))*normpdf(data_matrix(i,4)-(t1-t2),0,sqrt(v))) (line 36)
* normpdf(data_matrix(i, 2) – (t1 + t2), 0, sqrt(v)) …
^
Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), …
^^^^^^^^^^^^^^^^^^^^^^^
Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);
^^^^^^
Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, …
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integralCalc (line 113)
[q,errbnd] = vadapt(@minusInfToInfInvTransform,interval, …
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( …
^^^^^^^^^^^^^^^^^
Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( …
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);
^^^^^^
Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, …
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integralCalc (line 113)
[q,errbnd] = vadapt(@minusInfToInfInvTransform,interval, …
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integral2Calc (line 7)
[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Structural_Estimation>neg_log_lik_real (line 43)
pdf_value = integral2(f, -Inf, Inf, -Inf, Inf, ‘AbsTol’, 1e-6, ‘RelTol’, 1e-6);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Structural_Estimation>@(params)neg_log_lik_real(params,data_matrix) (line 58)
neg_log_lik_handle = @(params) neg_log_lik_real(params, data_matrix);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fminunc (line 233)
f = feval(funfcn{3},x,varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Structural_Estimation (line 62)
[estimated_params, fval, exitflag, output] = fminunc(neg_log_lik_handle, initial_params, options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Caused by:
Failure in initial objective function evaluation. FMINUNC cannot continue.

Related documentationDear community,
I am struggling with an error that my code produces. I am trying to do a maximum-likelihood estimation (MLE). I am importing CSV data; each row is an observation and an observation is four-dimensional with column names "bad_bought", "good_bought", "bad_notbought", "good_notbought".