Integrate function that calls a function

18 views (last 30 days)
Nils
Nils on 12 Aug 2014
Commented: Mike Hosea on 13 Aug 2014
I'm having trouble integrating a function that contains another function. To give an example:
value = X*func(y)
prob = lognpdf(y, muy, stdy)
int_func = @(y)(value(y)*prob(y))
i.e. I have a function value that depends on a variable y, which is log-normally distributed. The problem is that the multiplication in the first line does not work when done inside the integral function:
>> integral(int_func, 0, 5)
Error using _*_: Innder matrix dimensions must agree
As func(y) creates a matrix, the size of which does not line up with the size of X when it is not passed a scalar as input. It appears that integral() internally passes vectors to the function that is being integrated, which then causes the error. Is there any way around this?

Accepted Answer

Mike Hosea
Mike Hosea on 13 Aug 2014
Edited: Mike Hosea on 13 Aug 2014
Except when computing array-valued integrals, the INTEGRAL function requires that the integrand function be coded to accept arrays as input and return an array of the same size as output. We call this "vectorization". It is critical to performance in MATLAB. One rarely needs matrix multiplication * when defining an integrand function, rather element-wise multiplication .* .
If I understand you correctly, this is one of those cases where an intermediate quantity is an array, and yet the final result of evaluating the integrand function with a scalar is a scalar. We just need to understand that INTEGRAL is going to pass you an array, and it expects you to evaluate your integrand function at each element of the array. Given that, and the simplicity of the problem definition, we can go one of three ways with this.
One is to use 'ArrayValued',true. This is for calculating the integral of array-valued functions, usually vector-valued functions. Of course a scalar is a vector of length 1. It's the same computation, but in this mode INTEGRAND just passes scalars to the integrand function because it is expecting vectors (or any other fixed shape) back. So using 'ArrayValued',true with a scalar-valued problem essentially has the effect of turning off "vectorization" inside of the integral function. As a result, it will be slower than if we fixed the integrand to accept and return arrays, but the results should be valid.
Two is to use arrayfun. In this approach, you define your integrand f(y) however is natural for you, but you don't integrate f(y) directly, rather you pass this
fv = @(y)arrayfun(f,y);
to INTEGRAL. That is the syntax when f is a function_handle. If, on the other hand, you had an m-file named myfun.m, you would write
fv = @(y)arrayfun(@myfun,y);
These ARRAYFUN definitions mean "Take this function (the first input) and evaluate it element-by-element on entries of the following array (the second input).". It's like writing a loop. This is usually faster than 'ArrayValued',true on scalar-valued integrals.
Finally, in this case you might be able to write func in such a way that it would accept a row (or column) vector input and return an array of size 128-by-numel(y). Then X*func(y) would return a row-vector, which we would then reshape back to size y, i.e.
value = @(y)reshape(X*func(y(:).'),size(y))
prob = @(y)lognpdf(y, muy, stdy)
int_func = @(y)(value(y).*prob(y))
Note the use of elementwise multiplication .* rather than * in the definition of int_fun.
  3 Comments
Mike Hosea
Mike Hosea on 13 Aug 2014
I did a little example and got 18x speedup by rewriting func to accept a vector input and return a 128-by-numel(y) array. I don't know what's typical, but I'm not surprised by that number. It's significant.

Sign in to comment.

More Answers (1)

Honglei Chen
Honglei Chen on 12 Aug 2014
Edited: Honglei Chen on 12 Aug 2014
What is the dimensions of value and prob? What are the operation you want to do?
If the function is vector valued, integral can deal with it by setting 'ArrayValued' to true
  1 Comment
Nils
Nils on 13 Aug 2014
Hi Honglei,
thanks for your answer. Using the ArrayValued argument does give me an answer, although I'm not entirely sure it's the correct one (I don't have any analytical expression for the function I'm integrating, so verifying is not straightforward). Can you give an intuitive explanation of what ArrayValued does in this case?
Before I answer your question about the dimensions of my objects, let me note that I have made a mistake in my original post. The second line should have read:
prob = @(y)(lognpdf(y, muY, stdY))
so that both value and prob are functions of y. The dimensions are as follows:
>> size(X)
1 128
>> size(func(y)) % with y scalar
128 1
>> size(value)
1 1
>> size prob(y) % with y scalar
1 1
so that in the case of scalar y I simply evaluate value at y, which gives a scalar, and multiply by the probability of y occuring, which is also scalar. However, when I integrate using integral, it seems that a vector is passed to the func function in value, so that the multiplication X*func(y) fails. This is as the size of X is fixed, but func(y) changes shape when a vector is passed as an argument rather than a scalar.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!