How to find a solution using vpasolve or maybe other?

Question

Mikie 2018-2-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/381851-how-to-find-a-solution-using-vpasolve-or-maybe-other

评论： Walter Roberson 2018-6-10

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 个评论
显示 12更早的评论隐藏 12更早的评论

Mikie 2018-2-11

在 MATLAB Online 中打开

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

Mikie 2018-2-14

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

Mikie 2018-2-28

在 MATLAB Online 中打开

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

Mikie 2018-3-5

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

请先登录，再回答此问题。