Computation of sum/product of triple integral using function handle
4 次查看(过去 30 天)
显示 更早的评论
Consider a set of non-negative real values
,
and
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788470/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788475/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788480/image.png)
First, I want to compute f as a function of θ vector:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788485/image.png)
However, I am mainly interested in calculating the natural logarithm of this function. Namely:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788490/image.png)
I coded this for the case of a single
as:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788495/image.png)
% 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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/788495/image.png)
3 个评论
采纳的回答
Walter Roberson
2021-11-28
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.
2 个评论
Walter Roberson
2021-11-28
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])
更多回答(1 个)
Matt J
2021-11-4
编辑:Matt J
2021-11-4
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Whos 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!