Entering edit mode
@antillenicolaslausannenrc-bas-337
Last seen 10.0 years ago
Hi!
During background subtraction (with RMA), it is assumed that
the the observed intensities (O) are the sum of a signal(S)
and a noise(N) :
O = S + N
with S follows an exponential ditribution (parameter = alpha) and
N follows a normal distribution (parameters = mu and sigma).
I've tried to assess the quality of these parameters estimates by
simulation. For that, I've chosen the values of the three parameters
and I've generated a large data set (N values):
dataset <- rexp(N,alpha) + rnorm(N,mu,sigma)
After that, i've tried to estimate the values of the parameters like
in RMA:
bg.parameters(dataset, 2^12)
and I've been surprised to discover that the estimates weren't good at
all
For example, with alpha=0.0005, mu=500, lambda=75 and N=50000, I've
obtained
the following estimates:
alpha = 0.0029
mu = 776.64
sigma = 251.14
This encouraged me to find another way to estimate the parameters.
I've
found
it, but I need your opinion:
First, I've computed the theoretical distribution of O (convolution of
S and
N).
I've found the following expression:
f(o) = alpha * exp(alpha*(mu + 0.5*alpha*sigma^2 - o)) *
phi((o - mu
- alpha*sigma^2)/(sigma))
where phi is the cumulative distribution function of the normal
distribution.
Then I've computed the corresponding log-likelihood function. The
computation of maximum
log-likelihood is possible but difficult. That's why I propose to
solve it
in a different way.
We have: E[O] = E[S+N] = E[S]+E[N] = 1/alpha + mu = (appr.)
mean(intensities)
and Var[O] = Var[S+N] = Var[S] + Var[N] = 1/alpha^2 +
sigma^2 =
(appr.) var(intensities)
where we suppose that S and N are independant and where intensities is
the
vector of intensities (perfect match).
Then I suggest to make the following substitutions in the log-
likelihood
function:
mu ----> mean(intensities) - 1/alpha
sigma^2 ----> var(intensities) - 1/alpha^2
At this point, we have to maximize the log-likelihood which contains
just
one unknown parameter : alpha
For the moment, I suggest to find the maximum empirically (compute the
function for a large range of alpha)
In the example, I've computed the function for 5000 values of alpha
randomly
selected between 0.000005 and 0.05.
I've obtained the following estimates:
alpha = 0.000498
mu = 501.09
sigma = 113.84
It seems much better than the previous estimates.
I'm impatiently waiting for your comments and criticals.
Thank you in advance!
Nicolas
