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 v.

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 1x2, 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 documentation

Answer 1

[Torsten]

"integral2" usually calls the function to be integrated with matrices of the same size for t1 and t2 to make the necessary function evaluation more efficient.

You will have to use "arrayfun" to evaluate f only with scalar inputs or define f in a function (not as a function handle) and loop over the elements in t1 and t2.

If possible, check what you get from integral2. Usually - since the mass of a pdf is concentrated in a small region of IR^2 - numerical integration with integral2 gives bad results. That's why e.g. an extra function "mvncdf" exists to integrate "mvnpdf" using a specialized algorithm (instead of using integral2).

mu1 = 56;
mu2 = 11;
v1 = 119;
v2 = 131;
rho = -0.36;
v = 419;
    
% Mean vector
mu_vec = [mu1, mu2];
    
% Covariance matrix
cov_matrix = [v1, rho * sqrt(v1 * v2); rho * sqrt(v1 * v2), v2];

% Define the integrand for 2-dimensional integration
f = @(t1, t2) ( ...
    mvnpdf([t1, t2], mu_vec, cov_matrix) ...
    * normpdf(1 - (t1 - t2), 0, sqrt(v)) ...
    * normpdf(2 - (t1 + t2), 0, sqrt(v)) ...
    * normpdf(3 - (t1 + t2), 0, sqrt(v)) ...
    * normpdf(4 - (t1 - t2), 0, sqrt(v)) ...
    );
 F = @(t1,t2)arrayfun(@(t1,t2)f(t1,t2),t1,t2)

F = function_handle with value:
    @(t1,t2)arrayfun(@(t1,t2)f(t1,t2),t1,t2)

 F([2 5],[8 7])

ans = 1×2

1.0e-15 * 0.0624 0.2516

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

integral2(F, -Inf, Inf, -Inf, Inf, 'AbsTol', 1e-6, 'RelTol', 1e-6)

ans = 6.9476e-11

f(2,8)

ans = 6.2442e-17

f(5,7)

ans = 2.5163e-16

f([2 5],[8 7])

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

Error in solution>@(t1,t2)(mvnpdf([t1,t2],mu_vec,cov_matrix)*normpdf(1-(t1-t2),0,sqrt(v))*normpdf(2-(t1+t2),0,sqrt(v))*normpdf(3-(t1+t2),0,sqrt(v))*normpdf(4-(t1-t2),0,sqrt(v))) (line 16)
mvnpdf([t1, t2], mu_vec, cov_matrix) ...