Problem in solving equation with Gamma distribution parameters as unknowns

I am trying to solve the following sets of equations using the solve command.

gaminv(0.95,((2*y)-1),(1/(2*z)))-gaminv(0.05,((2*y)-1),(1/(2*z)))=6.73;

gaminv(0.45,((2*y)-1),(1/(2*z)))=7.19;

sqrt((Ia*(2*z)^((2*y)-1))/gamma((2*y)-1))=1.219;

It is giving an error of "Unable to prove '2*y - 1 < Inf & 0 < 2*y - 1 & 0 < 1/(2*z)' literally. Use 'isAlways' to test the statement mathematically."

I have also used the assume condition with the following lines

assume(0 < (2*y - 1) & (2*y - 1) < Inf);

assumeAlso(0 < (1/(2*z)));

Still same error. Is there any other way to solve them?

Answer 1

[Torsten]

I'd try it numerically. Better initial values for y,z and la might give convergence.

fun = @(y,z,la)[gaminv(0.95,y,z)-gaminv(0.05,y,z)-6.73;gaminv(0.45,y,z)-7.19;la/z^y-1.219^2*gamma(y)];
F = @(x)fun(x(1),x(2),x(3));
y0 = 1; z0 = 1; la0 = 1 ;
x0 = [y0,z0,la0];
sol = fsolve(F,x0,optimset('MaxIter',2500000,'MaxFunEvals',2500000));

Solver stopped prematurely.

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

norm(F(sol))

ans = 0.2273

y = 2*sol(1)-1

y = 25.2578

z = 1/(2*sol(2))

z = 0.8860

la = sol(3)

la = 5.3964e+05

[OUGATA MUKHERJEE]

Thank you for the quik reply. It worked with some minor modifications.

[Torsten]

To derive the "original" y parameter from your equations, I made a small mistake.

You have to use

y = (sol(1)+1)/2

instead of

y = 2*sol(1)-1

in the above code.