Integral2 seems to substitute non-scalar values of variable into integrand. Why?
21 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Lukas
am 13 Sep. 2024 um 20:35
Bearbeitet: Torsten
am 14 Sep. 2024 um 0:42
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
0 Kommentare
Akzeptierte Antwort
Torsten
am 13 Sep. 2024 um 22:25
Bearbeitet: Torsten
am 14 Sep. 2024 um 0:42
"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([2 5],[8 7])
integral2(F, -Inf, Inf, -Inf, Inf, 'AbsTol', 1e-6, 'RelTol', 1e-6)
f(2,8)
f(5,7)
f([2 5],[8 7])
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Matrix Computations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!