How to maximize this function?? Please help me.

Hi, I have some question about how to maximize this function

L = n*log(2/sigma) + sum(log(normpdf((W-mu_0)/sigma))) + sum(log(normcdf((W-mu_0)/sigma))) - (3*pi^2/32)*log(1+8*lambda^2/pi^2)

at (sigma,lambda). So, what I did is I took the derivative w.r.t sigma and lambda and then solve the those two partial derivative equations (w.r.t. sigma and lambda) = 0. I used fsolve and vpasolve but the solutions are quite different and some observations gave no solutions. Can anyone help me how to deal with this?

global W n m1 mu_0
N = 1000;
n = 5;
lambda = 0;
mu_0 = 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;
[sigmaRMM,lambdaRMM] = [1.0483,-5.3559];
%%%%%%%%%%%%%%%%%%%%%%%%%%
fun = @root2dmu0;

function F = root2dmu0(x)
global W n m1 mu_0
F(1) = (-1) - (x(2)/n)*sum((W - mu_0)/x(1))*normpdf(x(2)*(W-mu_0)/x(1))/normcdf(x(2)*(W-mu_0)/x(1));
F(2) = -2*(3*pi^2/32)*(8/pi^2)*x(2)/(1+8*x(2)^2/pi^2) + sum((W - mu_0)/x(1))*normpdf(x(2)*(W-mu_0)/x(1))/normcdf(x(2)*(W-mu_0)/x(1));

    x0 = [sigmaRMM,lambdaRMM]
    options = optimset('Display','off');
    Sol = fsolve(fun,x0,options);
    sigma = Sol(1)
    lambda = Sol(2)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    syms lambda_Hat sigma_Hat

    eqn_sigma = (-1) + (1/n).*sum(((W-mu_0).^2)/(sigma_Hat.^2)) - (lambda_Hat./(n*sigma_Hat)).*sum((W-mu_0).*normpdf((lambda_Hat.*(W-mu_0)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_0)./sigma_Hat)));
    eqn_lambda = ((-3*pi^2.*lambda_Hat)./(2.*(pi.^2+(8.*lambda_Hat.^2)))) + sum(((W-mu_0)./sigma_Hat).*normpdf((lambda_Hat.*(W-mu_0)./sigma_Hat))./normcdf((lambda_Hat.*(W-mu_0)./sigma_Hat)));

    S = vpasolve([eqn_sigma==0,eqn_lambda==0],[sigma_Hat,lambda_Hat],[sigmaRMM,lambdaRMM]);
    sigma = S.sigma_Hat
    lambda = S.lambda_Hat

[John D'Errico]

Please stop posting the same question. This question is identical to what you posted last time. The code you show here was in the comments to your last question.

[Walter Roberson]

That previous question does not appear to be available anymore.

Answer 1

[Walter Roberson]

You have

F(1) = (-1) - (x(2)/n)*sum((W - mu_0)/x(1))*normpdf(x(2)*(W-mu_0)/x(1))/normcdf(x(2)*(W-mu_0)/x(1));

but W is a 5 x 1 vector, so the denominator normcdf(x(2)*(W-mu_0)/x(1)) is a 5 x 1 vector and the subexpression in the numerator normpdf(x(2)*(W-mu_0)/x(1)) is a 5 x 1 vector. That gives you a 5 x 1 vector "/" a 5 x 1 vector. The / operator is matrix right division, where A/B is approximately the same as A * pinv(B) . (5x1) / (5x1) gives a 5 x 5 output, which cannot be assigned into the single location F(1) .

If you change the / to be ./ so that you do element-by-element division instead of matrix division, then you would have (5x1) ./ (5x1) giving a 5 x 1 output, which cannot be assigned to the single location F(1).

If you are trying to fit coefficients, finding the single coefficient that best explains the two 5x1 with respect to each other, then you would use the \ operator: 5 x 1 \ 5 x 1 would give 1 x 1.

... but I don't think that is what you are trying to do.

I suspect that you should be generating a single random value each time, instead of 5 values.

I also suspect that you are programming in Octave instead of MATLAB, as your code is not valid MATLAB code.