Computation of sum/product of triple integral using function handle

Consider a set of non-negative real values , and .

First, I want to compute f as a function of θ vector:

However, I am mainly interested in calculating the natural logarithm of this function. Namely:

I coded this for the case of a single as:

% Inputs:
t       = 15;
sigma   = 0.25;
y       = 1.0;

h  = @(x1,x2,x3)       x1.*(t - x2).^x3;
g  = @(x1,x2,x3)       exp(-0.5.*((y - h(x1,x2,x3))./sigma).^2 );        
f  = @(x1,x2,x3,theta) g(x1,x2,x3).*theta(1).*theta(2).*theta(3).*theta(4);            
J1 = @(theta)          integral3(@(x1,x2,x3) f(x1,x2,x3,theta), 0, 1, 0, t, 0.3, 1);
J2 = @(theta)          log(J1(theta));

...but I don't know how to do it for the general case of multiple . Any idea please?

[Matt J]

Note that the theta-dependence is trivial. Because the product over the thetas can be factored out of the triple integral, the function is really,

logf=n*sum(log(theta)) + constant

If this is a loglikelihood you are minimizing, the constant term is irrelevant. It is independent of theta and doesn't influence the minimization.

[DIMITRIS GEORGIADIS]

Hi Matt. First thank you for your answer. I think I understand what you are saying, but I am really interested in evaluating this constant. Could you please modify the code according to your understanding? (Maybe I have to use the arrayfun somewhere but I am not sure).

% Inputs:
t       = 15;
sigma   = 0.25;
y       = [1.2 1.0]';

h  = @(x1,x2,x3)       x1.*(t - x2).^x3;
g  = @(x1,x2,x3)       exp(-0.5.*((y - h(x1,x2,x3))./sigma).^2 );        
f  = @(x1,x2,x3,theta) g(x1,x2,x3).*theta(1).*theta(2).*theta(3).*theta(4);            
J1 = @(theta)          integral3(@(x1,x2,x3) f(x1,x2,x3,theta), 0, 1, 0, t, 0.3, 1);
J2 = @(theta)          sum(log(J1(theta)));
% Check:
q  = J2([0.062 1.21 7.5 0.65])

[DIMITRIS GEORGIADIS]

@Ameer Hamza @Walter Roberson can you help me on that?

Answer 1

[Walter Roberson]

h  = @(x1,x2,x3)       x1.*(t - x2).^x3;
g  = @(x1,x2,x3)       exp(-0.5.*((y - h(x1,x2,x3))./sigma).^2 );        
f  = @(x1,x2,x3,theta) g(x1,x2,x3).*theta(1).*theta(2).*theta(3).*theta(4);            
J1 = @(theta)          integral3(@(x1,x2,x3) f(x1,x2,x3,theta), 0, 1, 0, t, 0.3, 1);

When you use integral3(), the function must accept three arrays the same size, and must return an output the same size. It must perform element-wise multiplication.

So integral3 is going to pass arrays for x1, x2, x3, and @(x1,x2,x3) f(x1,x2,x3,theta) must return something the same size. The result of integral3() will be a scalar.

This has the direct implication that the size of the output of @(x1,x2,x3) f(x1,x2,x3,theta) must not depend upon the size of theta. And it also means that

J2 = @(theta)          sum(log(J1(theta)));

the J1 is going to return a scalar, and it does not make sense to take sum() of a scalar.

Your code is obviously hoping that you can use integral3() to return something the same size as theta, but integral3() cannot do that.

To return something the same size as theta, you are going to have to use arrayfun, such as

J2 = @(Theta) sum(log(arrayfun(J1, Theta)));

Or perhaps it is not theta you are concerned about...

Look again at the sequence. If y were non-scalar then g would not return something the same size as the x values, and that would mean that f would return something that was not the same size as the x values, and that would mean that @(x1,x2,x3) f(x1,x2,x3,theta) would return something that was not the same size as the x values, but returning something the same size is a requirement of integral3() . So y cannot be a non-scalar.

Perhaps you were hoping that integral3() would return one value for each y... but it cannot do that. You would have to arrayfun() over the Y values and pass the current y into the functions.

[Walter Roberson]

t       = 15;
sigma   = 0.25;
y       = [1.2 1.0]';
h  = @(x1,x2,x3)         x1.*(t - x2).^x3;
g  = @(x1,x2,x3,Y)       exp(-0.5.*((Y - h(x1,x2,x3))./sigma).^2 );        
f  = @(x1,x2,x3,theta,Y) g(x1,x2,x3,Y).*theta(1).*theta(2).*theta(3).*theta(4);            
J1 = @(theta,Y)          integral3(@(x1,x2,x3) f(x1,x2,x3,theta,Y), 0, 1, 0, t, 0.3, 1);
J2 = @(theta)            sum(log(arrayfun(@(Y) J1(theta,Y), y)));
% Check:
q  = J2([0.062 1.21 7.5 0.65])

q = -0.5948

[DIMITRIS GEORGIADIS]

Thank you very much @Walter Roberson. This is exactly what I was looking for!

Answer 2

[Matt J]

If you can do this for a single y_i, then it should be trivial for multiple y_i. Note that the y_i can be factored out of the integral as well, so that

logf=n*sum(log(theta)) -sum(y)/2 + constant

where the constant is obtained by doing the integral with all theta=1 and y=0. Since you say you can already do these, the problem is solved.