integral2 with arrayvalued

I have a fairly complex function that I want to do double integration as the following:

r1 = 0.0185;

r2 = 0.0172;

sig1 = 0.2222;

sig2 = 0.1450;

sigF = 0.1013;

rhot = -0.5112;

rho = 0.7908;

T = 2;

muX = (r1 - rhot*sig1*sigF - 0.5*(sig1^2))*T;

muY = (r2 - 0.5*(sig2^2))*T;

sigX = sqrt((sig1^2)*T);

sigY = sqrt((sig2^2)*T);

mu = [muX, muY]; 

sig = [sigX^2, rho*sigX*sigY; rho*sigX*sigY, sigY^2]; 

ymin = @(x) log(exp(x)+0.05);

fun = @(x,y) (exp(x) - exp(y))' * reshape(mvnpdf([x(:),y(:)],mu,sig), size(x));

D = integral2(fun, -Inf, Inf, ymin, Inf , 'RelTol',1e-10);

The error says that I have to set ArrayValued to true. The problems are 1. I don't know the size of the numerical data 2. There isn't a parameter to set in integral2. How can I solve this problem?

Answer 1

[Torsten]

function main
r1 = 0.0185;
r2 = 0.0172;
sig1 = 0.2222;
sig2 = 0.1450;
sigF = 0.1013;
rhot = -0.5112;
rho = 0.7908;
T = 2;
muX = (r1 - rhot*sig1*sigF - 0.5*(sig1^2))*T;
muY = (r2 - 0.5*(sig2^2))*T;
sigX = sqrt((sig1^2)*T);
sigY = sqrt((sig2^2)*T);
mu = [muX, muY]; 
sig = [sigX^2, rho*sigX*sigY; rho*sigX*sigY, sigY^2]; 
ymin = @(x) log(exp(x)+0.05);
D = integral2(@(x,y)fun(x,y,mu,sig), -Inf, Inf, ymin, Inf , 'RelTol',1e-10);

function f = fun(x,y,mu,sig)
for i=1:numel(x)
  vec = [x(i) y(i)];
  f(i) = (exp(x(i)-exp(y(i))*mvnpdf(vec,mu,sig);
end

Best wishes

Torsten.

[N/A]

I want to know a bit more about why my problem exist and how this change solved it.

[Torsten]

You will need pointwise multiplication instead of vector multiplication in your function "fun". I think Walters's answer shows the difference.

Best wishes

Torsten.

Answer 2

[Walter Roberson]

The function will be called with two arrays the same size, and must return a value the same size.

Suppose the array size is 3 x 2. Then you have

fun = @(x,y) (exp(x) - exp(y))' * reshape(mvnpdf([x(:),y(:)],mu,sig), size(x));

The (exp(x) - exp(y))' part would be 2 x 3 because of the transpose. The reshape(mvnpdf([x(:),y(:)],mu,sig), size(x)) part is going to be 3 x 2, the same size as the input.

You now have a 2 x 3 matrix multiplied by a 3 x 2. The result is going to be 2 x 2, but the size of result needed is 3 x 2.

In general if the input is m x n, your output is going to be n x n instead of m x n.

Perhaps you need

fun = @(x,y) (exp(x) - exp(y)) .* reshape(mvnpdf([x(:),y(:)],mu,sig), size(x));

[N/A]

Thank you!