Computation of sum/product of triple integral using function handle

3 ビュー (過去 30 日間)
DIMITRIS GEORGIADIS
DIMITRIS GEORGIADIS 2021 年 11 月 3 日
コメント済み: DIMITRIS GEORGIADIS 2021 年 11 月 29 日
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?
  3 件のコメント
DIMITRIS GEORGIADIS
DIMITRIS GEORGIADIS 2021 年 11 月 3 日
編集済み: DIMITRIS GEORGIADIS 2021 年 11 月 3 日
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
DIMITRIS GEORGIADIS 2021 年 11 月 4 日
編集済み: DIMITRIS GEORGIADIS 2021 年 11 月 4 日
@Ameer Hamza @Walter Roberson can you help me on that?

サインインしてコメントする。

採用された回答

Walter Roberson
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
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])
q = -0.5948
DIMITRIS GEORGIADIS
DIMITRIS GEORGIADIS 2021 年 11 月 29 日
Thank you very much @Walter Roberson. This is exactly what I was looking for!

サインインしてコメントする。

その他の回答 (1 件)

Matt J
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.

カテゴリ

Help Center および File ExchangeWhos についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by