How to do a double integral on a function with 3 variables, so the answer is dependent on the 3 variable

20 views (last 30 days)
Trond Oesten
Trond Oesten on 13 Mar 2015
Commented: Mike Hosea on 16 Mar 2015
Hi,
I have a function with 3 variables. I want to do a double integral on my function so it's dependent on the last variable. The variables I want to integrate is u and k so my function is dependent on x. I tryed using integral2, but that didn't work. Is it an easy way I can do this?
Thanks!
My script:
clc; clear all; close all;
C = [23.875 15.75281; 15.75281 93.9842];
mu_U = 1788.2058;
mu_K = 70.8489;
sigma_U = sqrt(C(1,1));
sigma_k = sqrt(C(2,2));
ubot = mu_U-4*sigma_U;
utop = mu_U+4*sigma_U;
kbot = mu_K-4*sigma_k;
ktop = mu_K+4*sigma_k;
mu = [mu_U;mu_K];
fun = @(x,u,k) wblpdf(x,u,k)*mvnpdf([u;k],mu,C);

Accepted Answer

Mike Hosea
Mike Hosea on 16 Mar 2015
Edited: Mike Hosea on 16 Mar 2015
Not sure what your problem with INTEGRAL2 was, but you're using matrix multiplication in your definition of FUN, which probably isn't right--I doubt you want to multiply a matrix of WBLPDF values by a matrix of MVNPDF values, and you're stacking u and k as different rows when you input to MVNPDF when I think what you should be doing is passing them as different columns. Assuming you meant to integrate over u and k first and to have a function of x coming out, this might be the sort of thing you want:
C = [23.875 15.75281; 15.75281 93.9842];
mu_U = 1788.2058;
mu_K = 70.8489;
sigma_U = sqrt(C(1,1));
sigma_k = sqrt(C(2,2));
ubot = mu_U-4*sigma_U;
utop = mu_U+4*sigma_U;
kbot = mu_K-4*sigma_k;
ktop = mu_K+4*sigma_k;
mu = [mu_U;mu_K];
fun = @(x,u,k) wblpdf(x,u,k).*reshape(mvnpdf([u(:),k(:)],mu(:).',C),size(u));
qfun_scalar = @(x)integral2(@(u,k)fun(x,u,k),ubot,utop,kbot,ktop,'AbsTol',1e-10,'RelTol',1e-6);
qfun = @(x)arrayfun(qfun_scalar,x);
I just threw the tolerances in there in case you wanted to play with them.
Now you can use qfun to evaluate for x values (scalar or arrays), or you can integrate over ranges of x, or whatever.
>> qfun(1750)
ans =
0.006956904221022
>> integral(qfun,1600,1800)
ans =
0.789977339847138
  2 Comments
Mike Hosea
Mike Hosea on 16 Mar 2015
Iqfun_scalar = @(a,b)integral(qfun,a,b);
As the name implies, that function is only set up for scalar a and scalar b. Now to vectorize it, I'll add a little trick using ONES to make it do scalar expansion.
Iqfun = @(a,b)arrayfun(Iqfun_scalar,a.*ones(size(b)),b.*ones(size(a)));
>> Iqfun([1600,1700],1800)
ans =
0.789977339847138 0.758915659832745
Of course you're doing a triple integral there, a couple of them, so it isn't blazingly fast.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!