이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How to find a solution using vpasolve or maybe other?
조회 수: 1 (최근 30일)
이전 댓글 표시
Hi, Can anyone please help me how to do this? I want to solve these 3 equations but MATLAB didn't give me any solutions sometimes when the data is simulated.
N = 100;
n=5;
lambda = 1;
for i = 1:N
U = normrnd(0,1,n,1);
V = normrnd(0,1,n,1);
W = (lambda/sqrt(1+lambda^2))*abs(U) + (1/sqrt(1+lambda^2))*V;
m1 = sum(W)/n;
m2 = sum((W-m1).^2)/n;
m3 = sum((W-m1).^3)/n;
A = (m2^3)/(m3^2);
if A > ((pi-2)^3)/(2*(4-pi)^2)
theta_Hat = (2/pi) + (A*(2/pi)*((4/pi)-1)^2)^(1/3);
lambda_Hatt = sqrt(1./(theta_Hat-1));
else
lambda_Hatt = 0;
end
lambdaMM = sign(m3).*abs(lambda_Hatt);
sigmaMM = sqrt(m2./(1-((2/pi).*(lambdaMM.^2./(1+lambdaMM.^2)))));
muMM = m1 - (sigmaMM.*sqrt(2./pi).*lambdaMM./sqrt(1+lambdaMM.^2));
eqn_mu = ((m1 - mu_Hat)./sigma_Hat) - (lambda_Hat./n).*sum(normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_sigma = (-n./sigma_Hat) + (1/sigma_Hat.^3).*sum((W-mu_Hat).^2) - (lambda_Hat./sigma_Hat.^2).*sum((W-mu_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_lambda = ((-3*pi^2.*lambda_Hat)./(2.*(pi.^2+(8.*lambda_Hat.^2)))) + sum(((W-mu_Hat)./sigma_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
equations = [eqn_mu,eqn_sigma,eqn_lambda];
vars = [mu_Hat,sigma_Hat,lambda_Hat];
initials = [muMM,sigmaMM,lambdaMM];
[muML,sigmaML,lambdaML] = vpasolve(equations==0, vars, initials);
end
Please help. Thank you so much.
댓글 수: 14
Mikie
2018년 2월 11일
Sorry.
syms theta lambda_Hat mu_Hat sigma_Hat
N = 1000;
n=5;
lambda = 1;
for i = 1:N
U = normrnd(0,1,n,1);
V = normrnd(0,1,n,1);
W = (lambda/sqrt(1+lambda^2))*abs(U) + (1/sqrt(1+lambda^2))*V;
m1 = sum(W)/n;
m2 = sum((W-m1).^2)/n;
m3 = sum((W-m1).^3)/n;
A = (m2^3)/(m3^2);
if A > ((pi-2)^3)/(2*(4-pi)^2)
theta_Hat = (2/pi) + (A*(2/pi)*((4/pi)-1)^2)^(1/3);
lambda_Hatt = sqrt(1./(theta_Hat-1));
else
lambda_Hatt = 0;
end
lambdaMM = sign(m3).*abs(lambda_Hatt);
sigmaMM = sqrt(m2./(1-((2/pi).*(lambdaMM.^2./(1+lambdaMM.^2)))));
muMM = m1 - (sigmaMM.*sqrt(2./pi).*lambdaMM./sqrt(1+lambdaMM.^2));
eqn_mu = ((m1 - mu_Hat)./sigma_Hat) - (lambda_Hat./n).*sum(normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_sigma = (-n./sigma_Hat) + (1/sigma_Hat.^3).*sum((W-mu_Hat).^2) - (lambda_Hat./sigma_Hat.^2).*sum((W-mu_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_lambda = ((-3*pi^2.*lambda_Hat)./(2.*(pi.^2+(8.*lambda_Hat.^2)))) + sum(((W-mu_Hat)./sigma_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
equations = [eqn_mu,eqn_sigma,eqn_lambda];
vars = [mu_Hat,sigma_Hat,lambda_Hat];
initials = [muMM,sigmaMM,lambdaMM];
[muML,sigmaML,lambdaML] = vpasolve(equations==0, vars, initials);
end
Walter Roberson
2018년 2월 11일
Not sure why you are overwriting your entire output variables each time through the loop?
Walter Roberson
2018년 2월 11일
I find that the values of muML and lambdaML jump around a lot.
I put in a breakpoint and tested the equations at the first point that vpasolve() failed to return a value. It turns at that in that case, there is a solution (at least in limit form) of sigma_Hat = +/- infinity and lambda_Hat = 0 and mu_Hat = anything. If the other equations have the same structure then this would explain why there is so much jumping around with large values of sigma_Hat and why the mu_Hat value is pretty much meaningless: you are getting whichever side of 0 happens to be tried first.
Mikie
2018년 2월 14일
I am a little bit confused and do not understand. Sorry, however, how can I tell MATLAB to keep only when it gives us solutions for muML, sigmaML, lambdaML and discard those replications which do not give solutions. and then maybe just keep track of the number of the successful replication which is used to find MSE_muMM,MSE_muML,Bias_muMM,Bias_muML later.
Thank you so much.
syms theta lambda_Hat mu_Hat sigma_Hat Lpmle
N = 10000;
n=5;
lambda = 1;
MSE_muMM = 0;
MSE_muML = 0;
Bias_muML = 0;
Bias_muMM = 0;
for i = 1:N
U = normrnd(0,1,n,1);
V = normrnd(0,1,n,1);
W = (lambda/sqrt(1+lambda^2))*abs(U) + (1/sqrt(1+lambda^2))*V;
m1 = sum(W)/n;
m2 = sum((W-m1).^2)/n;
m3 = sum((W-m1).^3)/n;
A = (m2^3)/(m3^2);
if A > ((pi-2)^3)/(2*(4-pi)^2)
theta_Hat = (2/pi) + (A*(2/pi)*((4/pi)-1)^2)^(1/3);
lambda_Hatt = sqrt(1./(theta_Hat-1));
else
lambda_Hatt = 0;
end
lambdaMM = sign(m3).*abs(lambda_Hatt);
sigmaMM = sqrt(m2./(1-((2/pi).*(lambdaMM.^2./(1+lambdaMM.^2)))));
muMM = m1 - (sigmaMM.*sqrt(2./pi).*lambdaMM./sqrt(1+lambdaMM.^2));
Lpmle = n*log(2/sigma_Hat) + sum(log(normpdf((W-mu_Hat)./sigma_Hat))) + sum(log(normcdf((lambda_Hat*(W-mu_Hat)./sigma_Hat)))) - 3*pi.^2*log(1+8*lambda_Hat.^2/pi^2)/32;
eqn_mu = diff(Lpmle,mu_Hat);
eqn_sigma = diff(Lpmle,sigma_Hat);
eqn_lambda = diff(Lpmle,lambda_Hat);
equations = [eqn_mu,eqn_sigma,eqn_lambda];
vars = [mu_Hat,sigma_Hat,lambda_Hat];
initials = [muMM,sigmaMM,lambdaMM];
[muML,sigmaML,lambdaML] = vpasolve(equations==0,vars,initials);
MSE_muMM = ((i-1)/i)*MSE_muMM + (1/i)*muMM.^2;
MSE_muML = ((i-1)/i)*MSE_muML + (1/i)*muML.^2;
Bias_muML = ((i-1)/i)*Bias_muML + (1/i)*muML;
Bias_muMM = ((i-1)/i)*Bias_muMM + (1/i)*muMM;
end
MSE_muMM
MSE_muML
Bias_muMM
Bias_muML
Walter Roberson
2018년 2월 14일
Your equations have theoretical solutions with one of the parameters being either -inf or +inf, and no theoretical based solver can choose between those two. When you add limited precision of calculation, you add a large number of solutions where that variable merely assumes large magnitude leading to division of the numerators by values large enough that the sub expressions become negligible. vpasolve does not have a preference for perusing a positive value instead of a negative value: it uses a Newton type method that gives it projected values based upon the local conditions near where it starts, and if round off error happens to favour one side of 0 over the other then it is happy to persue that side. Therefore for nearly identical starting equations, it might run off towards negative infinity or positive infinity so if you get 9e7 for one iteration and -9e7 for the next it doesn't mean oscillating, it just means that there are multiple solutions with one being chosen by random factors.
I would suggest that you need to think more about constraints for your values.
Mikie
2018년 2월 26일
I have improved a bit of my program and I have a question. Why does the answer from
S = vpasolve([eqn_mu==0,eqn_sigma==0,eqn_lambda==0],[mu_Hat,sigma_Hat,lambda_Hat],[muMM,sigmaMM,lambdaMM]);
and
S = vpasolve([eqn_mu==0,eqn_sigma==0,eqn_lambda==0],[mu_Hat,sigma_Hat,lambda_Hat]);
solving by giving an initial value and no giving any intial value come out different?
Thank you so much.
clear all
syms theta lambda_Hat mu_Hat sigma_Hat
N = 10000
n=5
lambda = .1
k=0;
MSE_muMM = 0;
MSE_muML = 0;
Bias_muML = 0;
Bias_muMM = 0;
for i = 1:N
U = normrnd(0,1,n,1);
V = normrnd(0,1,n,1);
W = (lambda/sqrt(1+lambda^2))*abs(U) + (1/sqrt(1+lambda^2))*V;
m1 = sum(W)/n;
m2 = sum((W-m1).^2)/n;
m3 = sum((W-m1).^3)/n;
A = (m2^3)/(m3^2);
if A > ((pi-2)^3)/(2*(4-pi)^2)
theta_Hat = (2/pi) + (A*(2/pi)*((4/pi)-1)^2)^(1/3);
lambda_Hatt = sqrt(1./(theta_Hat-1));
else
lambda_Hatt = 0;
end
lambdaMM = sign(m3).*abs(lambda_Hatt);
sigmaMM = sqrt(m2./(1-((2/pi).*(lambdaMM.^2./(1+lambdaMM.^2)))));
muMM = m1 - (sigmaMM.*sqrt(2./pi).*lambdaMM./sqrt(1+lambdaMM.^2))
eqn_mu = ((m1 - mu_Hat)./sigma_Hat) - (lambda_Hat./n).*sum(normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_sigma = (-1) + (1/n).*sum(((W-mu_Hat).^2)/(sigma_Hat.^2)) - (lambda_Hat./(n*sigma_Hat)).*sum((W-mu_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
eqn_lambda = ((-3*pi^2.*lambda_Hat)./(2.*(pi.^2+(8.*lambda_Hat.^2)))) + sum(((W-mu_Hat)./sigma_Hat).*normpdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_Hat)./sigma_Hat)));
S = vpasolve([eqn_mu==0,eqn_sigma==0,eqn_lambda==0],[mu_Hat,sigma_Hat,lambda_Hat],[muMM,sigmaMM,lambdaMM]);
muML = S.mu_Hat;
if abs(muML-muMM) < 100
MSE_muMM = ((i-1)/i)*MSE_muMM + (1/i)*muMM.^2;
MSE_muML = ((i-1)/i)*MSE_muML + (1/i)*muML.^2;
Bias_muML = ((i-1)/i)*Bias_muML + (1/i)*muML;
Bias_muMM = ((i-1)/i)*Bias_muMM + (1/i)*muMM;
else
i=i;
k=k+1;
end
end
Walter Roberson
2018년 2월 26일
"Why does [...] solving by giving an initial value and no giving any intial value come out different?"
vpasolve() is not a global minimizer. It uses a quasi-Newton approach, so it is prone to getting stuck in local minima.
vpasolve() is explicitly a solver that uses a limited number of significant digits, software floating point, so if two different inputs give an output that is the same to within those limited number of digits, then it cannot tell the two apart.
Your equations involve exp(-expression / sigma_Hat^2) . For finite expression and large enough sigma_Hat, the exp() is going to underflow to 0, and once that starts happening, the system effectively cannot tell apart different values of sigma_Hat
Mikie
2018년 2월 27일
Could you suggest me what solver I should use to find the maxima of mu_Hat?
Thank you.
Walter Roberson
2018년 2월 28일
My tests hint that the solution for the equations has mu_Hat = min(W) - 5*eps, with all mu_Hat greater than min(W)+17*eps leading to NaN in the calculations, and with a quite steep gradient near min(W) - 5*eps
The values of the other variables have much much less influence on the calculation.
The tolerances are so tight that it must be assumed that a minor rearrangement of the calculation could give a quite different calculation for any one position.
Mikie
2018년 2월 28일
Thank you so much.
Now, I am trying to maximize this function
Lpmle = @(p)n*log(2/p(2)) + sum(log(normpdf((W-p(1))./p(2)))) + sum(log(normcdf((p(3)*(W-p(1))./p(2))))) - 3*pi.^2*log(1+8*p(3).^2/pi^2)/32;
but I can solve it. Could you please give me any suggestion how to do it? Thank you.
clear all
N = 10000;
n=5;
lambda = .1;
k=0;
MSE_muMM = 0;
MSE_muML = 0;
Bias_muML = 0;
Bias_muMM = 0;
for i = 1:N
U = normrnd(0,1,n,1);
V = normrnd(0,1,n,1);
W = (lambda/sqrt(1+lambda^2))*abs(U) + (1/sqrt(1+lambda^2))*V;
m1 = sum(W)/n;
m2 = sum((W-m1).^2)/n;
m3 = sum((W-m1).^3)/n;
A = (m2^3)/(m3^2);
if A > ((pi-2)^3)/(2*(4-pi)^2)
theta_Hat = (2/pi) + (A*(2/pi)*((4/pi)-1)^2)^(1/3);
lambda_Hatt = sqrt(1./(theta_Hat-1));
else
lambda_Hatt = 0;
end
lambdaMM = sign(m3).*abs(lambda_Hatt);
sigmaMM = sqrt(m2./(1-((2/pi).*(lambdaMM.^2./(1+lambdaMM.^2)))));
muMM = m1 - (sigmaMM.*sqrt(2./pi).*lambdaMM./sqrt(1+lambdaMM.^2));
p0 = [muMM,sigmaMM,lambdaMM];
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
Lpmle = @(p)n*log(2/p(2)) + sum(log(normpdf((W-p(1))./p(2)))) + sum(log(normcdf((p(3)*(W-p(1))./p(2))))) - 3*pi.^2*log(1+8*p(3).^2/pi^2)/32;
sol= fminunc(Lpmle,p0,options);
muML = sol(1);
if abs(muML-muMM) < 100
MSE_muMM = ((i-1)/i)*MSE_muMM + (1/i)*muMM.^2;
MSE_muML = ((i-1)/i)*MSE_muML + (1/i)*muML.^2;
Bias_muML = ((i-1)/i)*Bias_muML + (1/i)*muML;
Bias_muMM = ((i-1)/i)*Bias_muMM + (1/i)*muMM;
else
i=i;
k=k+1;
end
end
N
n
lambda
k
MSE_muMM
MSE_muML
Bias_muMM
Bias_muML
Walter Roberson
2018년 2월 28일
You are using fminunc for Lpmle, which is an unconstrained minimization. The Lpmle function perhaps has several minima, but one of them is p(2) -> +inf and p(3) = 0, which is a combination which gives -inf for Lpmle.
Mikie
2018년 3월 5일
Hi, Thank you so much. But can we tell MATLAB to ignore these following results? and go to the next replication to find the reasonable/correct answer.
fminunc stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 300 (the default value).
Or
fminunc stopped because it cannot decrease the objective function
along the current search direction.
Walter Roberson
2018년 6월 10일
No. You can test the exit flag and do something different if those cases are detected, but you cannot tell fminunc to detect the situations internally and try somewhere else.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Get Started with Symbolic Math Toolbox에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 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 (한국어)