Compute sensitivity coefficient from numerical gradient

조회 수: 6 (최근 30일)
mz123
mz123 2023년 11월 21일
편집: Matt J 2023년 11월 22일
Hello everybody. I have to rewrite a very basic R script in matlab. The script performs calculation of the sensitivity coefficients, i.e., the partial derivatives, for an example of law of propagation of uncertainty.
Here the R code:
# measurement function
mass2.x <- function(x) {
m.std = x[1]; dm.d = x[2];
diff = x[3]; dm.c = x[4]; dm.B = x[5]
m.std + dm.d + diff + dm.c + dm.B
}
require(numDeriv)
m.std = 10000.005; dm.d = 0.0; diff = mean(c(0.01,0.03,0.02))
dm.c = 0.0; dm.B = 0.0;
sens = grad(func=mass2.x,x=c(m.std,dm.d,diff,dm.c,dm.B))
m.std.u = 0.0225
dm.d.u = 0.015/sqrt(3); dm.c.u = 0.010/sqrt(3)
diff.u = 0.025/sqrt(3); dm.B.u = 0.010/sqrt(3)
m.x = mass2.x(c(m.std,dm.d,diff,dm.c,dm.B))
m.x.unc = sqrt(sum(sens^2*c(m.std.u,dm.d.u,diff.u,dm.c.u,dm.B.u)^2))
basically the measurement model is like
and from the example if I were to calculate the sensitivity coefficient I need to obtain
>> sens = [1 1 1 1 1]
since the model is linear in all the parameters.
Now I tried this approach but I'm not obtaining the same results for sens:
m_std = 10000.005;
dm_d = 0.0;
diff = mean([0.01,0.03,0.02]);
dm_c = 0.0;
dm_B = 0.0;
F = [m_std, dm_d, diff, dm_c, dm_B]
sens = gradient(F)
m_std_u = 0.0225;
dm_d_u = 0.015/sqrt(3);
dm_c_u = 0.010/sqrt(3);
diff_u = 0.025/sqrt(3);
dm_B_u = 0.010/sqrt(3);
m_x = sum(F)
m_x_unc = sqrt(sum(sens.^2 .* [m_std_u,dm_d_u,diff_u,dm_c_u,dm_B_u].^2))
Surely I'm missing something and I do not know if I'm using the function correctly.
Thank you for your help.

답변 (1개)

Matt J
Matt J 2023년 11월 22일
편집: Matt J 2023년 11월 22일
I have no background in R. However, if you are trying to take the numerical gradient of a function, there is no native Matlab function that will do that. However this FEX download will,
func=@(x) sum(x);
sens = gradest(func,[0,0,0,0,0])
sens = 1×5
1 1 1 1 1

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by