How to do a double integral on a function with 3 variables, so the answer is dependent on the 3 variable
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Trond Oesten
il 13 Mar 2015
Commentato: Mike Hosea
il 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);
0 Commenti
Risposta accettata
Mike Hosea
il 16 Mar 2015
Modificato: Mike Hosea
il 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 Commenti
Mike Hosea
il 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.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!