Parametric mixture model of survival time data in matlab

Good day every one, I am new to matlab and working on parametric mixture model of survival time data. I would like to fit a mixture of weibull and gamma distributions to estimate the mixing probabilities and the parameters of the distributions via Expectation Maximization Algorithm EM. If any one how has experience in modeling survival mixtures to help and put me through.

Answer 1

[Tom Lane]

Do you need to use the EM algorithm? The mle function may be able to handle a problem like this:

x = [gamrnd(2,2,50,1); wblrnd(20,2,100,1)];
f = @(x,p,a,b,c,d) p*gampdf(x,a,b) + (1-p)*wblpdf(x,c,d);
mle(x,'pdf',f,'start',[.5 2 2 10 2])

[nasrin]

Hi Tom,

I am using mle function to fit mixture of gamma distributions. so I do:

f = @(x,p,a,b,c,d) p*gampdf(x,a,b) + (1-p)*gampdf(x,c,d); mle(x,'pdf',f,'start',[.5 2 2 10 1])

but I get this error:

The PDF function returned negative or zero values

I think the problem is with choice of starting values. Is there any way to infer them from my data? Thank you :)

[okoth ochola]

Hi Tom, kindly explain the meaning of codes in words so that one can easily apply the same in other data sets, for example what does line one give, why have you used the commands and are the meaning of the values inside the brackets? What do they do. Kindly explain if you dont mind. Take for example gamrnd(2,2,50,1); what does it mean?

x = [gamrnd(2,2,50,1); wblrnd(20,2,100,1)];

f = @(x,p,a,b,c,d) p*gampdf(x,a,b) + (1-p)*wblpdf(x,c,d);

mle(x,'pdf',f,'start',[.5 2 2 10 2])

Answer 2

[Russ Adheaux]

I have a related question. I have a a 1D array, x: 0<x<5; 300 data points; sorted ascending. I want to fit Rayleigh distribution in the range 0<x<=xo and 2-parameter weibull in the range xo<x<5. Then I want to plot the CDF of these two distributions on a single (Rayleigh probability plot) plot. Tom's answer gave some clue, but I don't seem to get it working. Thanks.

[Tom Lane]

Maybe you could clarify some more. It sounds like you want to estimate a Rayleigh parameter and two Weibull parameters at least. Do you know xo or do you need to estimate that? I guess you want to truncate the Rayleigh to [0,xo] and the Weibull to [xo,5]. Do you also expect to estimate the amount of probability to allocate to the Rayleigh vs. the Weibull? Then is it true that you want a probability plot against an untruncated Rayleigh?

[Russ Adheaux]

Sorry for the missing info...

Yes, I know xo. (Once I get this working I will try to estimate it)

Yes, I want to truncate Rayleigh to [0,xo] and the 2-parameter Weibull to [xo,5].

Because I am specifying xo, I thought I am specifying probability allocation, in a way.

I would like to see 3 things on my Rayleigh Probability Plot 1)data 2) Straight line indicating Rayleigh distrib 3) Another bent line from xo onwards showing 2-parameter fitted Weibull distrib.

I expect my data to lie on Rayleigh line up to xo and then on the Weibull line from xo onwards.

[Tom Lane]

This is trickier. Here's an example that works sometimes. First generate samples from two truncated distributions:

a = raylrnd(1,100,1);
b = wblrnd(4,10,100,1);
x0 = 3;
x = [a(a<x0); b(b>x0 & b<5)];
[N,X] = hist(x);
bar(X,N/((X(2)-X(1))*length(x)),'hist')
shg

Define truncated Rayliegh and Weibull distributions and a mixture.

rpdf = @(x,p1) (x<x0) .* raylpdf(x,p1)/max(.001,raylcdf(x0,p1));
wpdf = @(x,p2,p3) (x>x0 & x<5) .* wblpdf(x,p2,p3) / max(.001,(wblcdf(5,p2,p3)-wblcdf(x0,p2,p3)));
mypdf = @(x,p1,p2,p3,p4) (exp(p4)/(1+exp(p4)))*rpdf(x,abs(p1)) + (1/(1+exp(p4)))*wpdf(x,abs(p2),abs(p3));

Use mle to fit the mixture.

opt = statset('mlecustom');
opt = statset(opt,'FunValCheck','off','MaxIter',1e6);
p = mle(x,'pdf',mypdf,'start',[1 4 10 0])

Superimpose on the data, along with the separate distributions for reference.

xx = linspace(0,5);
line(xx,rpdf(xx,p(1)),'Color','g')
line(xx,wpdf(xx,p(2),p(3)),'Color','c')
line(xx,mypdf(xx,p(1),p(2),p(3),p(4)),'Color','r')
legend('hist(x)','Rayleigh','Weibull','Combined')

Note that I parametrized the mixture using p4 rather than working with the mixture proportion directly, because p4 doesn't need to be constrained to the interval [0,1].

[Note that I edited an earlier version of this post where I neglected to restrict the Rayleigh and Weibull functions to their range.]