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.

回答(2 个)

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])

2 个评论

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 :)
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])

请先登录,再进行评论。

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.

3 个评论

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?
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.
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.]

请先登录,再进行评论。

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by