Matlab Numerical Integration - Syntax/Style Query

Hi,

Pretext: I have no formal background in matlab or maths in general, so apologies if any of the following doesn't make sense. I am also new to this forum, so apologies if this post is incorrect in any way. [n.b. I also posted this thread on physicsforums.com]

Context: Ok, so I wanted to find the response value for various parameter values, using the following equation:

http://www.nitric.net/%7Epete/pictures/pj.integral.GIF

[If you can't see the image, it is an integration with infinite bounds over an equation consisting of a Gaussian probability density function multiplied by the the square of a Gaussian cumulative density function plus one minus the square of another Gaussian cumulative density function]

I wasn't sure whether it was possible to do this using the matlab symbolic toolkit (syms), so I thought I'd take a crack at it using numerical integration, using the quad command. After much effort and confusion, I ended up with the following code.

Code:

---

mu=0;
sigma=1;
s=0;
alpha=0.4;
 
mu=sprintf('%f', mu);
sigma=sprintf('%f', sigma);
s=sprintf('%f', s);
alpha=sprintf('%f', alpha);
 
 
f=strcat('normpdf(x+',s,',',mu,',sqrt(2)*',sigma,').*(normcdf(x,',mu,',sqrt(2)*',alpha,'*',sigma,').^2 + (1 - normcdf(x,',mu,',sqrt(2)*',alpha,'*',sigma,')).^2)');

quad(f,-1000,1000)

---

This code *DOES* appear to work (i.e. yields the expected answers). But it strikes me as being more than a little crude. My questions therefore are as follows:

1. From a mathematical point of view, was using numerical integration with suitably wide bounds the right way to go about this problem?
2. From a programming point of view, is there a more elegant way to execute the above numerical integration? (i.e. a better approach than wrapping everything up as one big string which is then passed to the quad function)

Many thanks for your time,
Best,
pete

Quote:

---

Originally Posted by pete j

Hi,

Pretext: I have no formal background in matlab or maths in general, so apologies if any of the following doesn't make sense. I am also new to this forum, so apologies if this post is incorrect in any way. [n.b. I also posted this thread on physicsforums.com]

Context: Ok, so I wanted to find the response value for various parameter values, using the following equation:

http://www.nitric.net/%7Epete/pictures/pj.integral.GIF

[If you can't see the image, it is an integration with infinite bounds over an equation consisting of a Gaussian probability density function multiplied by the the square of a Gaussian cumulative density function plus one minus the square of another Gaussian cumulative density function]

I wasn't sure whether it was possible to do this using the matlab symbolic toolkit (syms), so I thought I'd take a crack at it using numerical integration, using the quad command. After much effort and confusion, I ended up with the following code.

Code:

---

mu=0;
sigma=1;
s=0;
alpha=0.4;
 
mu=sprintf('%f', mu);
sigma=sprintf('%f', sigma);
s=sprintf('%f', s);
alpha=sprintf('%f', alpha);
 
 
f=strcat('normpdf(x+',s,',',mu,',sqrt(2)*',sigma,').*(normcdf(x,',mu,',sqrt(2)*',alpha,'*',sigma,').^2 + (1 - normcdf(x,',mu,',sqrt(2)*',alpha,'*',sigma,')).^2)');

quad(f,-1000,1000)

---

This code *DOES* appear to work (i.e. yields the expected answers). But it strikes me as being more than a little crude. My questions therefore are as follows:

1. From a mathematical point of view, was using numerical integration with suitably wide bounds the right way to go about this problem?
2. From a programming point of view, is there a more elegant way to execute the above numerical integration? (i.e. a better approach than wrapping everything up as one big string which is then passed to the quad function)

Many thanks for your time,
Best,
pete

---

You should define the integrand as a function in a .m file or as an inline function definition and then pass that function to the quadrature routine (check how to handle the extra function arguments you will have)

The other thing you can do is transform variable of integration so that the range is finite.

CB

Update for the record:

In terms of Matlab code - Thanks, I have done as suggested (code below). Much less cluttered, thumbs up.

In terms of the mathematics - Shall have to plead ignorance here and say that I'm not sure how to go about transforming the variable of integration. However, the current (numerical) method appears to be adequate for both my current and forseeable purposes, so on this occasion I'm inclined to let sleeping dogs lie.

Thanks for taking the time,
Best,
pete



%e.g.
mu=0;
sigma=1;
s=0;
alpha=0.4;


%anon function handle method
dpc = @(x, mu, sigma, s, alpha)normpdf(x+s,mu,sqrt(2)*sigma).*(normcdf(x,mu ,sqrt(2)*alpha*sigma).^2 + (1 - normcdf(x,mu,sqrt(2)*alpha*sigma)).^2);
[ans,nfe]=quad(@(x)dpc(x,mu, sigma, s, alpha),-1000,1000)


%separate .m file method
[ans,nfe]=quad(@(x)dpc(x,mu, sigma, s, alpha),-1000,1000)
function y=dpc(x, mu, sigma, s, alpha)
y=normpdf(x+s,mu,sqrt(2)*sigma).*(normcdf(x,mu,sqr t(2)*alpha*sigma).^2 + (1 - normcdf(x,mu,sqrt(2)*alpha*sigma)).^2);
end