I don't understand how you want to define the hessian with respect to "par" for vectors x and x0.
Do you want to create 100 x 100 Hessians of size 2x2, one for each pair (x(i),x0(j)) ?
An issue with matlabFunction
5 次查看（过去 30 天）
显示 更早的评论
Hello,
I have a question and I think the best is to explain via an example. So, consider the following
mu=@(x,par)par(1).*x;sigma=@(x,par)par(2);
syms x0 x dt
par = sym('par', [1, 2]);
LL = -1 / 2 * (log(2 * pi * dt * sigma(x0, par)^2) + ((x - x0 - mu(x0, par) * dt) / (sigma(x0, par) * sqrt(dt)))^2);
h_LL = hessian(LL,par);
h_LL=matlabFunction(h_LL)
Now, assume that V is a 1-by-n given vector and x0=V(1:n-1), x=V(2:n), dt and par's are given scalars. for instance, assume that dt=1, x0=1:100, x=2:101, par1=par2=1. What I want to do is to calculate the following sum over the elements of V:
h_LL(1,1,1,1,2)+h_LL(1,1,1,2,3)+h_LL(1,1,1,3,4)+ ... + h_LL(1,1,1,100,101)
Of course I can calculate this sum using a for-loop. But, my actual 'n' and the number of my actual parameters are very big making this scheme impractical. The idea that comes to my mind is to find a way to remove the 'reshape' in the above function handle h_LL and supply the entire vevtors x and x0 rather than doing so element by element.
Unfortunately, I haver no control over matlabFunction and cannot tell it please give me a matrix rather than applying reshape to a scalar.
Any idea?
Thanks in advance,
Babak
采纳的回答
Dyuman Joshi
2023-12-11
编辑：Dyuman Joshi
2023-12-11
mu=@(x,par)par(1).*x;
sigma=@(x,par)par(2);
syms x0 x dt
par = sym('par', [1, 2]);
LL = -1 / 2 * (log(2 * pi * dt * sigma(x0, par)^2) + ((x - x0 - mu(x0, par) * dt) / (sigma(x0, par) * sqrt(dt)))^2);
temp = hessian(LL, par)
%Reshape the hessian as a column vector
h_LL = reshape(temp, [], 1)
%Convert to a function handle
fun=matlabFunction(h_LL);
%Sample inputs
DT=1; X0=1:100; X=2:101; PAR1=1; PAR2=1;
%Provide x0 and x as row vectors
out = fun(DT, PAR1, PAR2, X, X0)
%Convert the output to a 3D array,
%where each 2D matrix is a corresponding output
out = reshape(out, 2, 2, [])
%For comparison
subs(temp, [dt x0 x par], [DT X0(1) X(1) PAR1 PAR2])
subs(temp, [dt x0 x par], [DT X0(2) X(2) PAR1 PAR2])
2 个评论
更多回答（1 个）
Walter Roberson
2023-12-11
matlabFunction can produce vectorized code only if all of the following are true:
- the expression is scalar for scalar inputs
- the expression does not use piecewise()
- the expression does not involve an integral, unless the integral is one that has a closed-form solution. This applies regardless of whether the symbolic parameters appear in the expression to be integrated, or in the bounds of the integral (but for different reasons for those two situations)
- the expression does not use symmatrix (matlabFunction refuses to process those)
I might be forgetting some other cases.
The idea that comes to my mind is to find a way to remove the 'reshape' in the above function handle h_LL
You would have to do that by post processing, such as by having matlabFunction use the 'File' option and edit the file afterwards. Unless you wanted to risk putting a reshape.m on your MATLAB path, but that would be distinctly risky.
So what is the official solution for this kind of problem?
- You can loop; or
- You can use arrayfun() to apply the function to "corresponding" matrix inputs, which will automatically wrap the results; you would 'uniform', 0 to get a cell array of results. You could reshape() or permute() that cell to be 1 x 1 x n and then use cell2mat() -- or you could do cat(3, ThatCell{:}) to concatenate along the third dimension.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
您也可以从以下列表中选择网站：
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)