parameter estimation using censored data by fitting a Maximum likelihood to a differential equation

7 views (last 30 days)
SP
SP on 4 Apr 2019
Commented: Torsten on 4 Apr 2019
I have a population data (N) measured over certain time points(t). The rate of change of the population is modelled as a ODE as,
My intention is to estimate the parameter value of λ using the N(t) data.
However, the issue is, values below the detection limit of 10 are all set to 10 so that I am working with censored data.
Data is as follows:
t=[0 1 2 3 4 5 6];
N(t)=[25000 10 10 10 10 10 10];
I am trying to estimate the parameter value ofλ using these censored data by following the approach found in here
Since , converting the exponential decay to a exponential distribution funciton using a normalizing factor of c.
Then, , where =25,000.
So, I assumedand estimated the parameter using MLE. The code is as below:.
N=[25000,10,10,10,10,10,10];
cutOff=10;
c=(N<=cutOff);
start=mean(N);
opt = statset('FunValCheck','off');
[paramEsts,paramCIs] = mle(N, 'censoring',c,'pdf',@exppdf, 'cdf',@expcdf, ...
'start',start, 'lower',0,'options',opt)
figure
scatter(0:6,N)
hold on
lambda=paramEsts/N(1)
t=[0 1 2 3 4 5 6]
fittedVals=N(1)*(exp(-lambda.*t));
plot(0:6,fittedVals)
There are few thing that I don't understand about this approach.
  1. I am not sure if the approach that I ma taking to fit for censored data is correct at all?
  2. I am not sure if the way I converted the ODE solution to a exponential PDF is right?
  3. Since a normalizing factor is used, should the estimated value of lambda for the ODE system be what I get as paramEsts or is it paramEsts/N0?
  4. so, paramEsts value I get is = 2.506*10^4. So, is lambda=2.506*10^4 or lambda=2.506*10^4/25000?
  5. when I plot the data with the fitted values the fitted is above the censored data points. but shouldn’t it be the case that it should close or below the censored data? because censored data means, the equipment was not powerful enough to measure values below 10, so the actual value can be <=10 and not above it?
It would be a great help if someone can help me with these problems.

Answers (0)

Community Treasure Hunt

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

Start Hunting!