Code for the Maximum likelihood esitmate for the gamma distribution , both parameters unknown

In the book by [Klugmann: Loss Models] "HERE IS THE BOOK: BOOK" (pleasee see also the attachment Loss.jpg) [OR Here] on the page 383 I would like to see a command for Matlab that would compute the Maximum Likelihood Estimate

for the Gamma Distribution where both parameters α and θ are unknown , please see the attachment. I would like to see a command that would output for Example 15.4 these numbers: =2561.1 and =0.55616 for these data

[27,82,115,126,155,161,243,294,340,384,457,680,855,877,974,1193,1340,1884,2558,15743]

Answer 1

[Torsten]

format long
data = [27,82,115,126,155,161,243,294,340,384,457,680,855,877,974,1193,1340,1884,2558,15743];
p = mle(data,'Distribution','Gamma')

p = 1×2

1.0e+03 * 0.000556157797255 2.561143630512257

[Jan]

That's a very nice solution! Suppose that we haven't the command mle at our disposal, could be a sequence of commnads in Matlab then be given, using some sequence of basic numerical solution commands ?

Just for learning purposes and my curiosity.

[Torsten]

Define the log-likelihood function for the gamma distribution with the above data, differentiate it with respect to the parameters, set the two equations to 0 and solve them for the two optimal parameters.

[Jan]

Well, I have 2 supplemetray question to my OQ: how do I make 1 parameter of Gamma known: and the other should be computed by MLE principle. Also, how can I insert the censoring above 250 ?

[Jan]

@Torsten I'm a complete beginner, I cannot differentiate nor solve for roots. Could you kindly give me a MWE ?

[Torsten]

how do I make 1 parameter of Gamma known: and the other should be computed by MLE principle.

Not possible as far as I know. You will have to program this on your own following the way I described ( differentiating the maximum likelihood function with respect to the unknown parameter, setting the result to 0 and solving for the unknown parameter).

Also, how can I insert the censoring above 250 ?

What do you mean ? Maybe

data(data>250) = []

?

[Jan]

Perhaps my last query on this example: From where can I derive this number -162.29 , what should I take for x for the computed and ?

[Jan]

@Torsten Well, could you write for me at least this: Differentiate this function

x+sin x

and solve it for x=0 ?

Now I cannot differentiate and solve it for roots! :-(

>> solve(x-1,x*y+2)

ans =

struct with fields:

x: [1×1 sym]

y: [1×1 sym]

>>

[Torsten]

syms alpha theta x
data = sym('data',[1 20]);
f = 1/gamma(alpha) * 1/(theta^alpha) *x^(alpha-1) * exp(-x/theta);
fp = prod(subs(f,x,data))

fp = 

logfp = log(fp)

logfp = 

dlogfpdalpha = diff(logfp,alpha)

dlogfpdalpha = 

dlogfpdtheta = diff(logfp,theta)

dlogfpdtheta = 

Data = [27,82,115,126,155,161,243,294,340,384,457,680,855,877,974,1193,1340,1884,2558,15743];
dlogfpdalpha = subs(dlogfpdalpha,data,Data)

dlogfpdalpha = 

dlogfpdtheta = subs(dlogfpdtheta,data,Data)

dlogfpdtheta = 

sol = vpasolve([dlogfpdalpha==0,dlogfpdtheta==0],[alpha,theta],[0.5 2500])

sol = struct with fields:
    alpha: 0.55615779737149188594827727939715
    theta: 2561.143629976936064991373664248

[Jan]

Thank you for the detailed answer.

I wanted to impose that α will be between 0.5 and 0.9 while θ between 200 and 3000; this however

has failed:

sol = vpasolve([dlogfpdalpha==0,dlogfpdtheta==0],[alpha],[0.5 0.9],[theta],[200,3000])

Also this

sol = vpasolve([dlogfpdalpha==0,dlogfpdtheta==0],[alpha,theta],[0.5 1500])

gives theta more than 1500, namely 2561.14

[Torsten]

sol = vpasolve([dlogfpdalpha==0,dlogfpdtheta==0],[alpha,theta],[0.5 1500])

gives theta more than 1500, namely 2561.14

Read in the documentation of "vpasolve" what the third argument given to the solver means.

Hint: It doesn't impose any constraint on the solution.

And you cannot impose any constraints on the solution.

You have a system of two equations in two unknowns. This is solved by

alpha: 0.55615779737149188594827727939715 
theta: 2561.143629976936064991373664248

No way to influence these values in general.