
現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
Minimization problem with integral constraint
6 ビュー (過去 30 日間)
古いコメントを表示
Hello, I'm working with a 2D numerical density profile.
. I have a set of radius
a maximum radius R and I want to find the best fit to the data r. I know that I can use maximum likelihood or another method, but I have problems with the constraints for
, because I require that




At first I tried with bins and adjusted the curve with cftool, but I need more precision. So I want to use minimization with that constraint.
Thank you so much.
採用された回答
Matt J
2022 年 5 月 9 日
編集済み: Matt J
2022 年 5 月 9 日
Perhaps you could reparametrize the curve as,

which automatically satisfies the constraint for any b and c. Moreoever, since this form has only two unknown parameters, it should be relatively easy to do a parameter sweep to find at least a good initial guess of b and c.
7 件のコメント
Esteban Garcia
2022 年 5 月 9 日
編集済み: Esteban Garcia
2022 年 5 月 9 日
Hello Matt, your answer was extremely helpful!.. I already did the sweep and found c=-1.295 and b=0.76 and it makes sense.
Is there a method to get the minimization faster?, because I need to apply it to a lot of galaxy clusters. I've seen and tried some methods in the documentation but none of them seemed to accept variables as exponents..
For example, for the first galaxy I minimized this function for a and b
F=-(b^(706*c)*(1/b^2 + 1)^c*(4/b^2 + 1)^(2*c)*(1/(4*b^2) + 1)^(2*c)*(9/b^2 + 1)^(2*c)*(9/(4*b^2) + 1)^c*(1/(16*b^2) + 1)^c*(9/(16*b^2) + 1)^(5*c)*(25/(4*b^2) + 1)^c*(9/(25*b^2) + 1)^(4*c)*(16/(25*b^2) + 1)^(2*c)*(25/(16*b^2) + 1)^c*(49/(16*b^2) + 1)^c*(49/(25*b^2) + 1)^c*(64/(25*b^2) + 1)^(2*c)*(81/(25*b^2) + 1)^c*(9/(100*b^2) + 1)^c*(121/(16*b^2) + 1)^c*(121/(25*b^2) + 1)^(3*c)*(49/(100*b^2) + 1)^(3*c)*(169/(16*b^2) + 1)^(2*c)*(169/(25*b^2) + 1)^c*(121/(100*b^2) + 1)^c*(169/(100*b^2) + 1)^(2*c)*(324/(25*b^2) + 1)^c*(9/(400*b^2) + 1)^c*(49/(400*b^2) + 1)^c*(361/(100*b^2) + 1)^(5*c)*(81/(400*b^2) + 1)^(2*c)*(121/(400*b^2) + 1)^(4*c)*(169/(400*b^2) + 1)^(3*c)*(4/(625*b^2) + 1)^c*(529/(100*b^2) + 1)^(3*c)*(9/(625*b^2) + 1)^(2*c)*(16/(625*b^2) + 1)^c*(36/(625*b^2) + 1)^c*(49/(625*b^2) + 1)^c*(64/(625*b^2) + 1)^(2*c)*(289/(400*b^2) + 1)^(3*c)*(81/(625*b^2) + 1)^c*(121/(625*b^2) + 1)^c*(361/(400*b^2) + 1)^c*(144/(625*b^2) + 1)^(3*c)*(196/(625*b^2) + 1)^(2*c)*(729/(100*b^2) + 1)^c*(441/(400*b^2) + 1)^c*(256/(625*b^2) + 1)^(2*c)*(289/(625*b^2) + 1)^(2*c)*(529/(400*b^2) + 1)^(2*c)*(324/(625*b^2) + 1)^(2*c)*(361/(625*b^2) + 1)^c*(961/(100*b^2) + 1)^c*(441/(625*b^2) + 1)^(4*c)*(484/(625*b^2) + 1)^c*(529/(625*b^2) + 1)^(3*c)*(729/(625*b^2) + 1)^c*(961/(400*b^2) + 1)^(2*c)*(784/(625*b^2) + 1)^(2*c)*(841/(625*b^2) + 1)^(2*c)*(1089/(400*b^2) + 1)^(2*c)*(1024/(625*b^2) + 1)^(3*c)*(1369/(400*b^2) + 1)^c*(1156/(625*b^2) + 1)^(2*c)*(1296/(625*b^2) + 1)^(4*c)*(1521/(400*b^2) + 1)^c*(1521/(625*b^2) + 1)^c*(1849/(400*b^2) + 1)^c*(1764/(625*b^2) + 1)^(2*c)*(1849/(625*b^2) + 1)^c*(1/(2500*b^2) + 1)^c*(9/(2500*b^2) + 1)^c*(49/(2500*b^2) + 1)^(2*c)*(1936/(625*b^2) + 1)^c*(121/(2500*b^2) + 1)^c*(169/(2500*b^2) + 1)^c*(2116/(625*b^2) + 1)^(3*c)*(289/(2500*b^2) + 1)^(2*c)*(2401/(400*b^2) + 1)^c*(361/(2500*b^2) + 1)^c*(2304/(625*b^2) + 1)^c*(441/(2500*b^2) + 1)^c*(2601/(400*b^2) + 1)^c*(2401/(625*b^2) + 1)^(3*c)*(529/(2500*b^2) + 1)^c*(2601/(625*b^2) + 1)^c*(729/(2500*b^2) + 1)^(4*c)*(2704/(625*b^2) + 1)^(2*c)*(841/(2500*b^2) + 1)^(2*c)*(961/(2500*b^2) + 1)^c*(2916/(625*b^2) + 1)^(3*c)*(1089/(2500*b^2) + 1)^(4*c)*(3249/(400*b^2) + 1)^c*(1369/(2500*b^2) + 1)^c*(3249/(625*b^2) + 1)^(2*c)*(1521/(2500*b^2) + 1)^(3*c)*(3481/(625*b^2) + 1)^(3*c)*(1849/(2500*b^2) + 1)^c*(3844/(625*b^2) + 1)^c*(4096/(625*b^2) + 1)^c*(4356/(625*b^2) + 1)^c*(2601/(2500*b^2) + 1)^c*(4489/(625*b^2) + 1)^c*(4624/(625*b^2) + 1)^c*(5041/(625*b^2) + 1)^c*(3249/(2500*b^2) + 1)^c*(5184/(625*b^2) + 1)^c*(3481/(2500*b^2) + 1)^c*(3721/(2500*b^2) + 1)^c*(4489/(2500*b^2) + 1)^c*(5041/(2500*b^2) + 1)^(2*c)*(5329/(2500*b^2) + 1)^c*(7396/(625*b^2) + 1)^(2*c)*(5929/(2500*b^2) + 1)^c*(1/(10000*b^2) + 1)^c*(9/(10000*b^2) + 1)^(2*c)*(7569/(2500*b^2) + 1)^c*(81/(10000*b^2) + 1)^(2*c)*(289/(10000*b^2) + 1)^(2*c)*(361/(10000*b^2) + 1)^c*(529/(10000*b^2) + 1)^c*(961/(10000*b^2) + 1)^(3*c)*(1089/(10000*b^2) + 1)^(2*c)*(8649/(2500*b^2) + 1)^(2*c)*(1369/(10000*b^2) + 1)^c*(1521/(10000*b^2) + 1)^(3*c)*(1681/(10000*b^2) + 1)^c*(1849/(10000*b^2) + 1)^c*(9409/(2500*b^2) + 1)^(3*c)*(2209/(10000*b^2) + 1)^(4*c)*(9801/(2500*b^2) + 1)^c*(2401/(10000*b^2) + 1)^(2*c)*(2601/(10000*b^2) + 1)^(2*c)*(10201/(2500*b^2) + 1)^c*(2809/(10000*b^2) + 1)^(5*c)*(10609/(2500*b^2) + 1)^(2*c)*(3481/(10000*b^2) + 1)^c*(3721/(10000*b^2) + 1)^c*(3969/(10000*b^2) + 1)^(4*c)*(11881/(2500*b^2) + 1)^(6*c)*(4489/(10000*b^2) + 1)^c*(4761/(10000*b^2) + 1)^c*(12321/(2500*b^2) + 1)^(2*c)*(5041/(10000*b^2) + 1)^c*(12769/(2500*b^2) + 1)^(2*c)*(5329/(10000*b^2) + 1)^c*(13689/(2500*b^2) + 1)^c*(6561/(10000*b^2) + 1)^c*(6889/(10000*b^2) + 1)^c*(14641/(2500*b^2) + 1)^c*(7569/(10000*b^2) + 1)^(3*c)*(7921/(10000*b^2) + 1)^(3*c)*(8649/(10000*b^2) + 1)^c*(16641/(2500*b^2) + 1)^(2*c)*(9409/(10000*b^2) + 1)^c*(17161/(2500*b^2) + 1)^c*(17689/(2500*b^2) + 1)^c*(10201/(10000*b^2) + 1)^(2*c)*(10609/(10000*b^2) + 1)^(2*c)*(11449/(10000*b^2) + 1)^c*(11881/(10000*b^2) + 1)^c*(12321/(10000*b^2) + 1)^c*(12769/(10000*b^2) + 1)^c*(20449/(2500*b^2) + 1)^c*(13689/(10000*b^2) + 1)^c*(14161/(10000*b^2) + 1)^(3*c)*(14641/(10000*b^2) + 1)^(2*c)*(22801/(2500*b^2) + 1)^c*(16129/(10000*b^2) + 1)^c*(16641/(10000*b^2) + 1)^(2*c)*(17161/(10000*b^2) + 1)^(3*c)*(17689/(10000*b^2) + 1)^c*(19321/(10000*b^2) + 1)^c*(19881/(10000*b^2) + 1)^c*(20449/(10000*b^2) + 1)^(2*c)*(28561/(2500*b^2) + 1)^c*(21609/(10000*b^2) + 1)^(2*c)*(23409/(10000*b^2) + 1)^c*(24649/(10000*b^2) + 1)^c*(25281/(10000*b^2) + 1)^(3*c)*(25921/(10000*b^2) + 1)^c*(26569/(10000*b^2) + 1)^c*(27889/(10000*b^2) + 1)^(2*c)*(28561/(10000*b^2) + 1)^c*(29241/(10000*b^2) + 1)^c*(31329/(10000*b^2) + 1)^c*(32041/(10000*b^2) + 1)^c*(33489/(10000*b^2) + 1)^c*(35721/(10000*b^2) + 1)^c*(36481/(10000*b^2) + 1)^c*(37249/(10000*b^2) + 1)^c*(39601/(10000*b^2) + 1)^c*(40401/(10000*b^2) + 1)^(2*c)*(41209/(10000*b^2) + 1)^c*(42849/(10000*b^2) + 1)^(2*c)*(43681/(10000*b^2) + 1)^c*(44521/(10000*b^2) + 1)^c*(49729/(10000*b^2) + 1)^c*(51529/(10000*b^2) + 1)^(3*c)*(56169/(10000*b^2) + 1)^c*(57121/(10000*b^2) + 1)^c*(58081/(10000*b^2) + 1)^c*(59049/(10000*b^2) + 1)^(2*c)*(61009/(10000*b^2) + 1)^c*(63001/(10000*b^2) + 1)^(2*c)*(66049/(10000*b^2) + 1)^c*(69169/(10000*b^2) + 1)^c*(72361/(10000*b^2) + 1)^c*(73441/(10000*b^2) + 1)^c*(77841/(10000*b^2) + 1)^c*(80089/(10000*b^2) + 1)^(3*c)*(85849/(10000*b^2) + 1)^c*(95481/(10000*b^2) + 1)^c*(114921/(10000*b^2) + 1)^c*(116281/(10000*b^2) + 1)^c*(124609/(10000*b^2) + 1)^c*(137641/(10000*b^2) + 1)^c*(c + 1)^353)/(pi^353*((b^2 + 968562431735785/70368744177664)^(c + 1) - b^(2*c + 2))^353)
Again, thank you!
Matt J
2022 年 5 月 9 日
編集済み: Matt J
2022 年 5 月 9 日
How fast the minimization goes depends on how efficiently you have coded the objective function evaluations. You appear to be repeating many of the same intermediate computations in your one-line implementation. E.g. b^2 and 100*b^2 appear in multiple places.
Esteban Garcia
2022 年 5 月 9 日
編集済み: Esteban Garcia
2022 年 5 月 9 日
I did this:
%F is defined like sym
A=0;
for b=0.1:0.01:2
for c=-3.005:0.01:-0.015
n=vpa(subs(F),10);
if n<A
A=n;
B=b;
C=c;
end
end
end
So, you would recommend me to save b^2,100*b^2, etc. into variables and plug them in F?.
Matt J
2022 年 5 月 10 日
I would recommend implementing it as written in my original answer. No need to have the Symbolic Math Toolbox involved.
Esteban Garcia
2022 年 5 月 10 日
編集済み: Esteban Garcia
2022 年 5 月 10 日
Hello Matt, thank you for answer
I see, but evaluating directly didn't give me a result. When I evaluated the function into b=0.1:0.01:2, c=-3.005:0.01:-0.015
-(b^(706*c)*(1/b^2 + 1)^c*(4/b^2 + 1)^(2*c)*(1/(4*b^2) + 1)^(2*c)*(9/b^2 + 1)^(2*c)*(9/(4*b^2) + 1)^c*(1/(16*b^2) + 1)^c*(9/(16*b^2) + 1)^(5*c)*(25/(4*b^2) + 1)^c*(9/(25*b^2) + 1)^(4*c)*(16/(25*b^2) + 1)^(2*c)*(25/(16*b^2) + 1)^c*(49/(16*b^2) + 1)^c*(49/(25*b^2) + 1)^c*(64/(25*b^2) + 1)^(2*c)*(81/(25*b^2) + 1)^c*(9/(100*b^2) + 1)^c*(121/(16*b^2) + 1)^c*(121/(25*b^2) + 1)^(3*c)*(49/(100*b^2) + 1)^(3*c)*(169/(16*b^2) + 1)^(2*c)*(169/(25*b^2) + 1)^c*(121/(100*b^2) + 1)^c*(169/(100*b^2) + 1)^(2*c)*(324/(25*b^2) + 1)^c*(9/(400*b^2) + 1)^c*(49/(400*b^2) + 1)^c*(361/(100*b^2) + 1)^(5*c)*(81/(400*b^2) + 1)^(2*c)*(121/(400*b^2) + 1)^(4*c)*(169/(400*b^2) + 1)^(3*c)*(4/(625*b^2) + 1)^c*(529/(100*b^2) + 1)^(3*c)*(9/(625*b^2) + 1)^(2*c)*(16/(625*b^2) + 1)^c*(36/(625*b^2) + 1)^c*(49/(625*b^2) + 1)^c*(64/(625*b^2) + 1)^(2*c)*(289/(400*b^2) + 1)^(3*c)*(81/(625*b^2) + 1)^c*(121/(625*b^2) + 1)^c*(361/(400*b^2) + 1)^c*(144/(625*b^2) + 1)^(3*c)*(196/(625*b^2) + 1)^(2*c)*(729/(100*b^2) + 1)^c*(441/(400*b^2) + 1)^c*(256/(625*b^2) + 1)^(2*c)*(289/(625*b^2) + 1)^(2*c)*(529/(400*b^2) + 1)^(2*c)*(324/(625*b^2) + 1)^(2*c)*(361/(625*b^2) + 1)^c*(961/(100*b^2) + 1)^c*(441/(625*b^2) + 1)^(4*c)*(484/(625*b^2) + 1)^c*(529/(625*b^2) + 1)^(3*c)*(729/(625*b^2) + 1)^c*(961/(400*b^2) + 1)^(2*c)*(784/(625*b^2) + 1)^(2*c)*(841/(625*b^2) + 1)^(2*c)*(1089/(400*b^2) + 1)^(2*c)*(1024/(625*b^2) + 1)^(3*c)*(1369/(400*b^2) + 1)^c*(1156/(625*b^2) + 1)^(2*c)*(1296/(625*b^2) + 1)^(4*c)*(1521/(400*b^2) + 1)^c*(1521/(625*b^2) + 1)^c*(1849/(400*b^2) + 1)^c*(1764/(625*b^2) + 1)^(2*c)*(1849/(625*b^2) + 1)^c*(1/(2500*b^2) + 1)^c*(9/(2500*b^2) + 1)^c*(49/(2500*b^2) + 1)^(2*c)*(1936/(625*b^2) + 1)^c*(121/(2500*b^2) + 1)^c*(169/(2500*b^2) + 1)^c*(2116/(625*b^2) + 1)^(3*c)*(289/(2500*b^2) + 1)^(2*c)*(2401/(400*b^2) + 1)^c*(361/(2500*b^2) + 1)^c*(2304/(625*b^2) + 1)^c*(441/(2500*b^2) + 1)^c*(2601/(400*b^2) + 1)^c*(2401/(625*b^2) + 1)^(3*c)*(529/(2500*b^2) + 1)^c*(2601/(625*b^2) + 1)^c*(729/(2500*b^2) + 1)^(4*c)*(2704/(625*b^2) + 1)^(2*c)*(841/(2500*b^2) + 1)^(2*c)*(961/(2500*b^2) + 1)^c*(2916/(625*b^2) + 1)^(3*c)*(1089/(2500*b^2) + 1)^(4*c)*(3249/(400*b^2) + 1)^c*(1369/(2500*b^2) + 1)^c*(3249/(625*b^2) + 1)^(2*c)*(1521/(2500*b^2) + 1)^(3*c)*(3481/(625*b^2) + 1)^(3*c)*(1849/(2500*b^2) + 1)^c*(3844/(625*b^2) + 1)^c*(4096/(625*b^2) + 1)^c*(4356/(625*b^2) + 1)^c*(2601/(2500*b^2) + 1)^c*(4489/(625*b^2) + 1)^c*(4624/(625*b^2) + 1)^c*(5041/(625*b^2) + 1)^c*(3249/(2500*b^2) + 1)^c*(5184/(625*b^2) + 1)^c*(3481/(2500*b^2) + 1)^c*(3721/(2500*b^2) + 1)^c*(4489/(2500*b^2) + 1)^c*(5041/(2500*b^2) + 1)^(2*c)*(5329/(2500*b^2) + 1)^c*(7396/(625*b^2) + 1)^(2*c)*(5929/(2500*b^2) + 1)^c*(1/(10000*b^2) + 1)^c*(9/(10000*b^2) + 1)^(2*c)*(7569/(2500*b^2) + 1)^c*(81/(10000*b^2) + 1)^(2*c)*(289/(10000*b^2) + 1)^(2*c)*(361/(10000*b^2) + 1)^c*(529/(10000*b^2) + 1)^c*(961/(10000*b^2) + 1)^(3*c)*(1089/(10000*b^2) + 1)^(2*c)*(8649/(2500*b^2) + 1)^(2*c)*(1369/(10000*b^2) + 1)^c*(1521/(10000*b^2) + 1)^(3*c)*(1681/(10000*b^2) + 1)^c*(1849/(10000*b^2) + 1)^c*(9409/(2500*b^2) + 1)^(3*c)*(2209/(10000*b^2) + 1)^(4*c)*(9801/(2500*b^2) + 1)^c*(2401/(10000*b^2) + 1)^(2*c)*(2601/(10000*b^2) + 1)^(2*c)*(10201/(2500*b^2) + 1)^c*(2809/(10000*b^2) + 1)^(5*c)*(10609/(2500*b^2) + 1)^(2*c)*(3481/(10000*b^2) + 1)^c*(3721/(10000*b^2) + 1)^c*(3969/(10000*b^2) + 1)^(4*c)*(11881/(2500*b^2) + 1)^(6*c)*(4489/(10000*b^2) + 1)^c*(4761/(10000*b^2) + 1)^c*(12321/(2500*b^2) + 1)^(2*c)*(5041/(10000*b^2) + 1)^c*(12769/(2500*b^2) + 1)^(2*c)*(5329/(10000*b^2) + 1)^c*(13689/(2500*b^2) + 1)^c*(6561/(10000*b^2) + 1)^c*(6889/(10000*b^2) + 1)^c*(14641/(2500*b^2) + 1)^c*(7569/(10000*b^2) + 1)^(3*c)*(7921/(10000*b^2) + 1)^(3*c)*(8649/(10000*b^2) + 1)^c*(16641/(2500*b^2) + 1)^(2*c)*(9409/(10000*b^2) + 1)^c*(17161/(2500*b^2) + 1)^c*(17689/(2500*b^2) + 1)^c*(10201/(10000*b^2) + 1)^(2*c)*(10609/(10000*b^2) + 1)^(2*c)*(11449/(10000*b^2) + 1)^c*(11881/(10000*b^2) + 1)^c*(12321/(10000*b^2) + 1)^c*(12769/(10000*b^2) + 1)^c*(20449/(2500*b^2) + 1)^c*(13689/(10000*b^2) + 1)^c*(14161/(10000*b^2) + 1)^(3*c)*(14641/(10000*b^2) + 1)^(2*c)*(22801/(2500*b^2) + 1)^c*(16129/(10000*b^2) + 1)^c*(16641/(10000*b^2) + 1)^(2*c)*(17161/(10000*b^2) + 1)^(3*c)*(17689/(10000*b^2) + 1)^c*(19321/(10000*b^2) + 1)^c*(19881/(10000*b^2) + 1)^c*(20449/(10000*b^2) + 1)^(2*c)*(28561/(2500*b^2) + 1)^c*(21609/(10000*b^2) + 1)^(2*c)*(23409/(10000*b^2) + 1)^c*(24649/(10000*b^2) + 1)^c*(25281/(10000*b^2) + 1)^(3*c)*(25921/(10000*b^2) + 1)^c*(26569/(10000*b^2) + 1)^c*(27889/(10000*b^2) + 1)^(2*c)*(28561/(10000*b^2) + 1)^c*(29241/(10000*b^2) + 1)^c*(31329/(10000*b^2) + 1)^c*(32041/(10000*b^2) + 1)^c*(33489/(10000*b^2) + 1)^c*(35721/(10000*b^2) + 1)^c*(36481/(10000*b^2) + 1)^c*(37249/(10000*b^2) + 1)^c*(39601/(10000*b^2) + 1)^c*(40401/(10000*b^2) + 1)^(2*c)*(41209/(10000*b^2) + 1)^c*(42849/(10000*b^2) + 1)^(2*c)*(43681/(10000*b^2) + 1)^c*(44521/(10000*b^2) + 1)^c*(49729/(10000*b^2) + 1)^c*(51529/(10000*b^2) + 1)^(3*c)*(56169/(10000*b^2) + 1)^c*(57121/(10000*b^2) + 1)^c*(58081/(10000*b^2) + 1)^c*(59049/(10000*b^2) + 1)^(2*c)*(61009/(10000*b^2) + 1)^c*(63001/(10000*b^2) + 1)^(2*c)*(66049/(10000*b^2) + 1)^c*(69169/(10000*b^2) + 1)^c*(72361/(10000*b^2) + 1)^c*(73441/(10000*b^2) + 1)^c*(77841/(10000*b^2) + 1)^c*(80089/(10000*b^2) + 1)^(3*c)*(85849/(10000*b^2) + 1)^c*(95481/(10000*b^2) + 1)^c*(114921/(10000*b^2) + 1)^c*(116281/(10000*b^2) + 1)^c*(124609/(10000*b^2) + 1)^c*(137641/(10000*b^2) + 1)^c*(c + 1)^353)/(pi^353*((b^2 + 968562431735785/70368744177664)^(c + 1) - b^(2*c + 2))^353)
or
vpa(-(b^(706*c)*(1/b^2 + 1)^c*(4/b^2 + 1)^(2*c)*(1/(4*b^2) + 1)^(2*c)*(9/b^2 + 1)^(2*c)*(9/(4*b^2) + 1)^c*(1/(16*b^2) + 1)^c*(9/(16*b^2) + 1)^(5*c)*(25/(4*b^2) + 1)^c*(9/(25*b^2) + 1)^(4*c)*(16/(25*b^2) + 1)^(2*c)*(25/(16*b^2) + 1)^c*(49/(16*b^2) + 1)^c*(49/(25*b^2) + 1)^c*(64/(25*b^2) + 1)^(2*c)*(81/(25*b^2) + 1)^c*(9/(100*b^2) + 1)^c*(121/(16*b^2) + 1)^c*(121/(25*b^2) + 1)^(3*c)*(49/(100*b^2) + 1)^(3*c)*(169/(16*b^2) + 1)^(2*c)*(169/(25*b^2) + 1)^c*(121/(100*b^2) + 1)^c*(169/(100*b^2) + 1)^(2*c)*(324/(25*b^2) + 1)^c*(9/(400*b^2) + 1)^c*(49/(400*b^2) + 1)^c*(361/(100*b^2) + 1)^(5*c)*(81/(400*b^2) + 1)^(2*c)*(121/(400*b^2) + 1)^(4*c)*(169/(400*b^2) + 1)^(3*c)*(4/(625*b^2) + 1)^c*(529/(100*b^2) + 1)^(3*c)*(9/(625*b^2) + 1)^(2*c)*(16/(625*b^2) + 1)^c*(36/(625*b^2) + 1)^c*(49/(625*b^2) + 1)^c*(64/(625*b^2) + 1)^(2*c)*(289/(400*b^2) + 1)^(3*c)*(81/(625*b^2) + 1)^c*(121/(625*b^2) + 1)^c*(361/(400*b^2) + 1)^c*(144/(625*b^2) + 1)^(3*c)*(196/(625*b^2) + 1)^(2*c)*(729/(100*b^2) + 1)^c*(441/(400*b^2) + 1)^c*(256/(625*b^2) + 1)^(2*c)*(289/(625*b^2) + 1)^(2*c)*(529/(400*b^2) + 1)^(2*c)*(324/(625*b^2) + 1)^(2*c)*(361/(625*b^2) + 1)^c*(961/(100*b^2) + 1)^c*(441/(625*b^2) + 1)^(4*c)*(484/(625*b^2) + 1)^c*(529/(625*b^2) + 1)^(3*c)*(729/(625*b^2) + 1)^c*(961/(400*b^2) + 1)^(2*c)*(784/(625*b^2) + 1)^(2*c)*(841/(625*b^2) + 1)^(2*c)*(1089/(400*b^2) + 1)^(2*c)*(1024/(625*b^2) + 1)^(3*c)*(1369/(400*b^2) + 1)^c*(1156/(625*b^2) + 1)^(2*c)*(1296/(625*b^2) + 1)^(4*c)*(1521/(400*b^2) + 1)^c*(1521/(625*b^2) + 1)^c*(1849/(400*b^2) + 1)^c*(1764/(625*b^2) + 1)^(2*c)*(1849/(625*b^2) + 1)^c*(1/(2500*b^2) + 1)^c*(9/(2500*b^2) + 1)^c*(49/(2500*b^2) + 1)^(2*c)*(1936/(625*b^2) + 1)^c*(121/(2500*b^2) + 1)^c*(169/(2500*b^2) + 1)^c*(2116/(625*b^2) + 1)^(3*c)*(289/(2500*b^2) + 1)^(2*c)*(2401/(400*b^2) + 1)^c*(361/(2500*b^2) + 1)^c*(2304/(625*b^2) + 1)^c*(441/(2500*b^2) + 1)^c*(2601/(400*b^2) + 1)^c*(2401/(625*b^2) + 1)^(3*c)*(529/(2500*b^2) + 1)^c*(2601/(625*b^2) + 1)^c*(729/(2500*b^2) + 1)^(4*c)*(2704/(625*b^2) + 1)^(2*c)*(841/(2500*b^2) + 1)^(2*c)*(961/(2500*b^2) + 1)^c*(2916/(625*b^2) + 1)^(3*c)*(1089/(2500*b^2) + 1)^(4*c)*(3249/(400*b^2) + 1)^c*(1369/(2500*b^2) + 1)^c*(3249/(625*b^2) + 1)^(2*c)*(1521/(2500*b^2) + 1)^(3*c)*(3481/(625*b^2) + 1)^(3*c)*(1849/(2500*b^2) + 1)^c*(3844/(625*b^2) + 1)^c*(4096/(625*b^2) + 1)^c*(4356/(625*b^2) + 1)^c*(2601/(2500*b^2) + 1)^c*(4489/(625*b^2) + 1)^c*(4624/(625*b^2) + 1)^c*(5041/(625*b^2) + 1)^c*(3249/(2500*b^2) + 1)^c*(5184/(625*b^2) + 1)^c*(3481/(2500*b^2) + 1)^c*(3721/(2500*b^2) + 1)^c*(4489/(2500*b^2) + 1)^c*(5041/(2500*b^2) + 1)^(2*c)*(5329/(2500*b^2) + 1)^c*(7396/(625*b^2) + 1)^(2*c)*(5929/(2500*b^2) + 1)^c*(1/(10000*b^2) + 1)^c*(9/(10000*b^2) + 1)^(2*c)*(7569/(2500*b^2) + 1)^c*(81/(10000*b^2) + 1)^(2*c)*(289/(10000*b^2) + 1)^(2*c)*(361/(10000*b^2) + 1)^c*(529/(10000*b^2) + 1)^c*(961/(10000*b^2) + 1)^(3*c)*(1089/(10000*b^2) + 1)^(2*c)*(8649/(2500*b^2) + 1)^(2*c)*(1369/(10000*b^2) + 1)^c*(1521/(10000*b^2) + 1)^(3*c)*(1681/(10000*b^2) + 1)^c*(1849/(10000*b^2) + 1)^c*(9409/(2500*b^2) + 1)^(3*c)*(2209/(10000*b^2) + 1)^(4*c)*(9801/(2500*b^2) + 1)^c*(2401/(10000*b^2) + 1)^(2*c)*(2601/(10000*b^2) + 1)^(2*c)*(10201/(2500*b^2) + 1)^c*(2809/(10000*b^2) + 1)^(5*c)*(10609/(2500*b^2) + 1)^(2*c)*(3481/(10000*b^2) + 1)^c*(3721/(10000*b^2) + 1)^c*(3969/(10000*b^2) + 1)^(4*c)*(11881/(2500*b^2) + 1)^(6*c)*(4489/(10000*b^2) + 1)^c*(4761/(10000*b^2) + 1)^c*(12321/(2500*b^2) + 1)^(2*c)*(5041/(10000*b^2) + 1)^c*(12769/(2500*b^2) + 1)^(2*c)*(5329/(10000*b^2) + 1)^c*(13689/(2500*b^2) + 1)^c*(6561/(10000*b^2) + 1)^c*(6889/(10000*b^2) + 1)^c*(14641/(2500*b^2) + 1)^c*(7569/(10000*b^2) + 1)^(3*c)*(7921/(10000*b^2) + 1)^(3*c)*(8649/(10000*b^2) + 1)^c*(16641/(2500*b^2) + 1)^(2*c)*(9409/(10000*b^2) + 1)^c*(17161/(2500*b^2) + 1)^c*(17689/(2500*b^2) + 1)^c*(10201/(10000*b^2) + 1)^(2*c)*(10609/(10000*b^2) + 1)^(2*c)*(11449/(10000*b^2) + 1)^c*(11881/(10000*b^2) + 1)^c*(12321/(10000*b^2) + 1)^c*(12769/(10000*b^2) + 1)^c*(20449/(2500*b^2) + 1)^c*(13689/(10000*b^2) + 1)^c*(14161/(10000*b^2) + 1)^(3*c)*(14641/(10000*b^2) + 1)^(2*c)*(22801/(2500*b^2) + 1)^c*(16129/(10000*b^2) + 1)^c*(16641/(10000*b^2) + 1)^(2*c)*(17161/(10000*b^2) + 1)^(3*c)*(17689/(10000*b^2) + 1)^c*(19321/(10000*b^2) + 1)^c*(19881/(10000*b^2) + 1)^c*(20449/(10000*b^2) + 1)^(2*c)*(28561/(2500*b^2) + 1)^c*(21609/(10000*b^2) + 1)^(2*c)*(23409/(10000*b^2) + 1)^c*(24649/(10000*b^2) + 1)^c*(25281/(10000*b^2) + 1)^(3*c)*(25921/(10000*b^2) + 1)^c*(26569/(10000*b^2) + 1)^c*(27889/(10000*b^2) + 1)^(2*c)*(28561/(10000*b^2) + 1)^c*(29241/(10000*b^2) + 1)^c*(31329/(10000*b^2) + 1)^c*(32041/(10000*b^2) + 1)^c*(33489/(10000*b^2) + 1)^c*(35721/(10000*b^2) + 1)^c*(36481/(10000*b^2) + 1)^c*(37249/(10000*b^2) + 1)^c*(39601/(10000*b^2) + 1)^c*(40401/(10000*b^2) + 1)^(2*c)*(41209/(10000*b^2) + 1)^c*(42849/(10000*b^2) + 1)^(2*c)*(43681/(10000*b^2) + 1)^c*(44521/(10000*b^2) + 1)^c*(49729/(10000*b^2) + 1)^c*(51529/(10000*b^2) + 1)^(3*c)*(56169/(10000*b^2) + 1)^c*(57121/(10000*b^2) + 1)^c*(58081/(10000*b^2) + 1)^c*(59049/(10000*b^2) + 1)^(2*c)*(61009/(10000*b^2) + 1)^c*(63001/(10000*b^2) + 1)^(2*c)*(66049/(10000*b^2) + 1)^c*(69169/(10000*b^2) + 1)^c*(72361/(10000*b^2) + 1)^c*(73441/(10000*b^2) + 1)^c*(77841/(10000*b^2) + 1)^c*(80089/(10000*b^2) + 1)^(3*c)*(85849/(10000*b^2) + 1)^c*(95481/(10000*b^2) + 1)^c*(114921/(10000*b^2) + 1)^c*(116281/(10000*b^2) + 1)^c*(124609/(10000*b^2) + 1)^c*(137641/(10000*b^2) + 1)^c*(c + 1)^353)/(pi^353*((b^2 + 968562431735785/70368744177664)^(c + 1) - b^(2*c + 2))^353))
The result is always -inf...
But using the symbolyc function I got numbers like,
-2.966948035e-508
I think that is because the precision
Esteban Garcia
2022 年 5 月 11 日
Hello Matt it is much faster now without symbolyc and with log. Thank you
N=length(Rproj);
R=max(Rproj);
A=1000000
B=0
C=0
syms b c
for j=1:N
F(j)=log(2*Rproj(j)*(1+(Rproj(j)/b)^2)^c/(((b^2+R^2)^(c+1)-b^(2*c+2))/(b^(2*c)*(c+1))));
end
E=-sum(F);
fstr=string(E);
fstr=replace(fstr,'b','b(k)');
fstr=replace(fstr,'c','c(j)');
b=0.01:0.001:0.8
c=-1.405:0.001:-0.705
for k=1:length(b)
for j=1:length(c)
n=eval(fstr);
if n<A
A=n;
B=b(k);
C=c(j);
end
end
end
その他の回答 (1 件)
Mitch Lautigar
2022 年 5 月 9 日
My suggestion is to use a smaller step size for <a,b,c> if you know what they are. Typically when you are trying to fix the curve, the only thing you can do is try to add in more datapoints. If you can provide some code, I can provide more feedback.
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- 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 (한국어)