Finding roots of a two variable function

Basically I have x which is a column vector with 350 rows of random data. I need to model this as a Birnbaum Sauders distribution and estimate its parameters. I can do this directly using the mle function, but I wonder if its possible to use fsolve. Below is my function file.

function f = BirnbaumSaunders(alpha,lamb)

F = (alpha.^2)-((1/350)*sum((lamb.*x)+(1./(lamb.*x))-2));

G = (-175.*lamb)+((1./(2.*(alpha.^2)))*sum(((lamb.*lamb.*x))-(1./x))+sum(lamb./((lamb.*x)+1)));

end

Here, I have to find the value of alpha and lamb for which both my functions F and G are zero. I have tried using the fsolve with initial alpha and lamb values as 0.65 and 2.13 (I got these values after using the mle function). But I'm getting an error using fsolve.

fun = @BirnbaumSaunders;
x = fsolve(fun,[0.65,2.13])

This is giving me an error. Is there any other way to find the values of alpha and lamb when F and G are set to 0?

Answer 1

[Torsten]

x = rand(10,1);
fun = @(z)BirnbaumSaunders(z(1),z(2),x);
z = fsolve(fun,[0.65,2.13])

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 2.000000e+02.

z = 1×2

0.1961 -1.1427

norm(fun(z))

ans = 0.1632

function f = BirnbaumSaunders(alpha,lamb,x)
  F = (alpha.^2)-((1/350)*sum((lamb.*x)+(1./(lamb.*x))-2));
  G = (-175.*lamb)+((1./(2.*(alpha.^2)))*sum(((lamb.*lamb.*x))-(1./x))+sum(lamb./((lamb.*x)+1)));
  f = [F;G];
end

[Walter Roberson]

x = rand(10,1);
fun = @(z)BirnbaumSaunders(z(1),z(2),x);
opt = optimoptions('fsolve', 'MaxFunctionEvaluations', 1e4);
z = fsolve(fun,[0.65,2.13], opt)

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient,
but the vector of function values is not near zero as measured by the
value of the function tolerance.

z = 1×2

0.3520 -0.9806

norm(fun(z))

ans = 0.3046

function f = BirnbaumSaunders(alpha,lamb,x)
  F = (alpha.^2)-((1/350)*sum((lamb.*x)+(1./(lamb.*x))-2));
  G = (-175.*lamb)+((1./(2.*(alpha.^2)))*sum(((lamb.*lamb.*x))-(1./x))+sum(lamb./((lamb.*x)+1)));
  f = [F;G];
end

[Walter Roberson]

In my test, the coefficients for collect(F) become rather large, beyond 10^5000. Perhaps in practice the input x are significally less than 1 ??

[Chaitanya]

yes, so the x is basically a column vector and some sample data of x is shown below. It typically has data ranging between 0 and 1. can you try to do this now with the below x data?

0.260594
0.50998
0.609437
0.365165
0.949793
0.590385
0.902765

[Chaitanya]

I got a correct value using this data and your code!!! Thank you very much for your time!

[Torsten]

You want to estimate two parameters of a distribution given seven realizations ? That's ... courageous.

Answer 2

[Walter Roberson]

format long g
x = [
    0.260594
    0.50998
    0.609437
    0.365165
    0.949793
    0.590385
    0.902765];
fun = @(z)BirnbaumSaunders(z(1),z(2),sym(x));
syms z [1 2]
F = simplify(fun(z))

F = 

sol = vpasolve(F,[0.65;2.13])

sol = struct with fields:
    z1: [16×1 sym]
    z2: [16×1 sym]

vpa(subs(z, sol))

ans = 

subs(F, sol)

ans = 

fun = @(z)BirnbaumSaunders(z(1),z(2),x);
opt = optimoptions('fsolve', 'MaxFunctionEvaluations', 1e4);
Z = fsolve(fun,[0.06,2], opt)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Z = 1×2

0.0630815494828271 2.00076854993031

function f = BirnbaumSaunders(alpha,lamb,x)
  F = (alpha.^2)-((1/350)*sum((lamb.*x)+(1./(lamb.*x))-2));
  G = (-175.*lamb)+((1./(2.*(alpha.^2)))*sum(((lamb.*lamb.*x))-(1./x))+sum(lamb./((lamb.*x)+1)));
  f = [F;G];
end