Unbiased estimator of modulation depth?

[f91jsw]

I have the following problem:
I measure a sinusoidally varying signal at a number of phase points. I then fit a sine curve to the data points using least squares. The fitting function looks like:
f = a1 + a2 cos(phi) + a3 sin(phi)
I want to evaluate the modulation depth from the measurements given by
m = sqrt(a2*a2 + a3*a3)/a1
Now the problem is that mean(m) is a biased estimator of m, and this becomes very significant when the signal-to-noise ratio is low. This would seem like a pretty common type of problem but I'm not a stat guy. Can someone please point me in the right direction, literature, web site, anything that could help me out here?
Johannes

 

[EnumaElish]

Why is mean(m) biased? Does it have to do with m being a nonlinear function of the a's?

I guess that would be right, E[m] = E[sqrt(a₂² + a₃²)/a₁] is not equal to M = sqrt(A₂² + A₃²)/A₁ where uppercase denotes true parameter value.

 

[EnumaElish]

You could try "brute force" algebra to mould your expression into a shape that includes m as one of the parameters. That way you'd have a direct estimator.

E.g., if f = a₁ + a₂x + a₃z and the statistic you were looking to estimate was (a₁+a₂)², then your left hand side variable would be (1+1/x+1/x²)f², because the intercept term in that expression is a₁²+2a₁a₂+a₂².

 

[EnumaElish]

Although admittedly your problem is harder because it involves a square root and division.

 

[Hurkyl]

I can't point at literature, but maybe I can help figure this out. I'm not an electrical engineering expert, so you'll have to bear with me.


I think the scenario you're describing is the following:

The signal S(t) is being transmitted, where S(t) is given by:

S(t) = A_1 + A_2 \cos (\phi t) + A_3 \sin (\phi t)

where \phi is a known quantity. (?)


The signal is being transmitted through a noisy channel, and you sample the signal at times t_1, t_2, t_3, \cdots, yielding the samples:

y_k = S(t_k) + e_k

Where e_k, the error on the k-th sample, is gaussian noise.


Your goal is to estimate the quantity m = \sqrt{A_2^2 + A_3^2} / A_1 from the samples y_k.


Is that correct?

 

[f91jsw]

Your goal is to estimate the quantity m = \sqrt{A_2^2 + A_3^2} / A_1 from the samples y_k.
Is that correct?

Yes, correct, only that I don't measure S directly as a function of time, just phase. The measurement is based on a homodyne technique so I measure a DC signal at each discrete phase point (which I choose myself). It makes no difference for the problem though.

The way I do it now is, for an ensemble of measurements, to average each of the parameters a_1, a_2 and a_3 first. It works because it effectively increases the signal-to-noise ratio so any bias becomes negligible. It's a bit awkward though because for technical reasons the phase varies between the measurements so I need to make an accurate reference measurement to determine the phase variation in order to be able to average the a_i's.

Thanks for the input so far everyone!

Johannes

 

[Hurkyl]

I'm going to move this over to the Electrical Engineering section; maybe we can drum up some more knowledgeable people that way!

 

[EnumaElish]

The way I do it now is, for an ensemble of measurements, to average each of the parameters a_1, a_2 and a_3 first. It works because it effectively increases the signal-to-noise ratio so any bias becomes negligible.

Why do you need to average the a's? I am not an engineer, but why can't you just put all the observations you have collected into one big vector and let the regression do the averaging automatically (i.e. implicitly)?

 

[f91jsw]

Why do you need to average the a's? I am not an engineer, but why can't you just put all the observations you have collected into one big vector and let the regression do the averaging automatically (i.e. implicitly)?

That would accomplish the same thing, except that the varying phase complicates things. But still, it doesn't solve the fundamental problem of bias.

Johannes

 

[EnumaElish]

You may be able to run a simulation to determine the magnitude and the direction of the bias. Then come up with a correction formula.

 

[EnumaElish]

Another approach would be to linearize the m expression (e.g. using Fourier series); then manipulate the f function to estimate the linearized expression directly.

 

[EnumaElish]

I meant Taylor series, not Fourier.

 

[f91jsw]

You may be able to run a simulation to determine the magnitude and the direction of the bias. Then come up with a correction formula.

I've actually already done that with a short Matlab script. That's the engineering solution but I'm fundamentally interested in a proper statistical treatment of the problem. It seems like something that must pop up often: you have an unbiased maximum likelihood estimate of some parameters, and then use those in a non-linear formula which introduces bias. I was hoping to find a general approach to that.

Johannes

 

[EnumaElish]

Given m = G(a), you can always obtain an unbiased estimator of m's mean if you can derive (or simulate) the distribution of m from those of a.

Example: a and b are two normal variables, distributed with N(A,sigma_1) and N(B,sigma_2) respectively. You'd like to calculate the expected value of x = ab. Let X be the true mean of x. You can take 1,000 random pairs (a,b) from the two normal distributions and calculate the product of each pair. This will yield 1,000 observations of x, whose average is an unbiased estimator of X. (I am guessing that you have already tried this.)

Or you may be able to work out the analytical solution.

 

[f91jsw]

Given m = G(a), you can always obtain an unbiased estimator of m's mean if you can derive (or simulate) the distribution of m from those of a.

Maybe I misunderstand you but I think my problem is that I want an unbiased estimator of m, not the mean of m.

Example: assume I have a variable a distributed N(A,sigma), and m = sqrt(a). The the true m is now sqrt(A). If I have 1000 observations of a and take sqrt(a) for each of them and then the mean, mean(m) will not be = sqrt(A).

Johannes

 

[EnumaElish]

My post was poorly written.

m = sqrt(a) is a random variable distributed with F_m(u) = Prob(sqrt(a) < u) = Prob(a < u^2) = F_a(u^2). The first moment of F_m is the expected value of the estimator, m.

To empirically estimate the expected value of a random variable, one repeatedly measures observations of the variable and computes the arithmetic mean of the results. This estimates the true expected value in an unbiased manner ...