14
$\begingroup$

I am trying to extract a signal, $X$, from a noisy source of data, $Z$:

$Z = X + Y$

where I know that $Y\sim N(0,\sigma)$, and I know the standard deviation $\sigma$. I am starting with a surrogate distribution to the type of thing that I actually will be faced with. This is defined:

 fSurvivalGompertzDistRand[α_, β_] := ProbabilityDistribution[(1/((E^(α/β) Gamma[
      0, α/β])/β) E^(((1 - E^(t β)) α)/β)), {t, 0, ∞}]

If I draw independent samples from the above using:

data = RandomVariate[fSurvivalGompertzDistRand[0.016, 0.65], {20000}];

The histogram of results looks like:

enter image description here

I then add some gaussian noise to this signal, and plot a histogram of the results:

dataNoise = data + RandomVariate[NormalDistribution[0, 2.5], {20000}];
Histogram[dataNoise, 20]

Which results in a Histogram that looks like:

enter image description here

I would like to try to extract out the original signal, and have been looking into using a Gaussian filter. I have tried the following, but it doesn't seem to work:

Histogram[GaussianFilter[dataNoise, {1, 2.5}]]

In particular, I'm not sure how to choose the parameters of the filter, to best extract out the original signal. I presume that $\sigma$ in the Mathematica documentation is the $\sigma$ above? I have tried a large number of combinations of the parameters in the filter, but none seem any good at extracting out the original signal.

I have a fear that the amount of noise I am adding means that it is basically going to be impossible to reconstruct the original signal. The histogram of new data looks basically gaussian, so I fear the noise is just swamping the signal.

Does anyone have any ideas here? The thing I am most interested in extracting is the PDF of $X$, so if there are any methods that bypass the need to get $X$, then I'd also be interested in hearing them. Sorry, I am a complete novice in signal processing.

edit 1: Note, I am not interested in estimating the parameters of the distribution I define above - I am only using it as an example of the types of distribution that I might encounter in practice. I want to be able to estimate the PDF of $X$ for an arbitrary distribution.

$\endgroup$
  • $\begingroup$ Doesn't obtaining the estimates of the coefficients from your previous post (mathematica.stackexchange.com/questions/108052/…) get you the underlying signal? In other words you get estimates of $\alpha$ and $\beta$ which with the functional form of the Gompertz distribution gives the "signal". $\endgroup$ – JimB Feb 23 '16 at 23:23
  • $\begingroup$ Hi Jim, I can estimate these, although the estimates are very noisy. I was trying a completely different way of approaching the problem which is agnostic to the form of the distribution of $X$. This is because in practice, the age distribution will be unknown. Essentially, I want to try to directly infer its PDF - which from my previous question represents the age distribution. From this I can then calculate the hazard rate. Does that make sense? Not sure if I'm explaining myself well. Best, Ben $\endgroup$ – ben18785 Feb 23 '16 at 23:44
  • $\begingroup$ If the precision of the estimates from the method of moments is not sufficient, it is unlikely that maximum likelihood or any other frequentist method can do much better in this case. And even a Bayesian approach would require you to put some awfully precise and accurate priors on the three parameters to improve the precision. I could be wrong but you can't change the laws of statistics (or physics) to paraphrase a future engineer. In any event, this question might be better served on Cross Validated. $\endgroup$ – JimB Feb 24 '16 at 2:36
  • 1
    $\begingroup$ I am inclined to agree, however the method I am proposing here is sort of orthogonal to what you mention. I am trying to directly estimate the PDF, not its parameters. I fear that either way the method may be doomed due to the small size of the population age structure relative to the bandwidth of the gaussian noise. If I don't get a reasonable answer here, I will repost on Cross Validated, or somewhere else (CV never seems to work for me.) Best, Ben $\endgroup$ – ben18785 Feb 24 '16 at 2:39
  • $\begingroup$ @ben18785 Did you post this question also to Cross Validated? $\endgroup$ – Anton Antonov Jul 23 '16 at 18:06
19
$\begingroup$

The procedure below assumes that the original distribution $X$ (the "signal") is non-Gaussian, and $Y$ is Gaussian (normally distributed noise.)

General procedure

The procedure is as follows:

  1. Find a function $F$ that applied to a collection of real numbers produces one value (say, 0) for normally distributed data and other different values for non-normally distributed data.

  2. Pick a sample $\{s_i\}_{i=1}^{n}$ of the noised data with a relatively small number of $n$ points (say $n \in [20,40]$).

  3. Formulate an optimization problem with $n$ variables $\{v_i\}_{i=1}^{n}$ that maximizes $ \lvert F(s-v) \rvert$ subject to constrains that would hold if $\{v_i\}_{i=1}^{n}$ come from a normal distribution (with known parameters).

  4. Solve the optimization problem several times with different samples and accumulate the sets $\{s_i - v_i\}_{i=1}^{n}$ and $\{v_i\}_{i=1}^{n}$. Monitor the $\chi^2$ test over $\{v_i\}_{i=1}^{n}$.

  5. Reconstruct the original distribution (the "signal") CDF and PDF using quantiles of the union of $\{s_i - v_i\}_{i=1}^{n}$ from all optimization runs.

Step details

Non-Gaussianity measure

First, we are going to adopt exess kurtosis as a measure for non-Gaussianity. For this we are going to rewrite excess kurtosis as

ExKurtosis[inp_] := CentralMoment[inp, 4] - 3 CentralMoment[inp, 2]^2

For points comming from the Normal Distribution ExKurtosis is close to 0:

In[1139]:= ExKurtosis[NormalDistribution[a, b]]
Out[1139]= 0

So, excess kurtosis close to 0 means Gaussianity, excess kurtosis significantly larger than 0 means non-Gaussianity.

Other non-Gaussiany measures exist with better properties (theoretical justification, robustness, speed of computation). See this article "Independent Component Analysis: A Tutorial" .

Constraints for normally distributed noise

We should come up with constraints which would hold if the values given to the variables are normally distributed. Since we know the mean and the standard deviation of the noise we can write up several such constraints based on the properties of Normal Distribution. (Mean, StandardDeviation, Kurtosis, etc.)

Constraints from "signal" knowlege

We can add constraints coming up from our knowledge of the distribution that is noised.

From the examples in the question we can add the constraints $\{s_i - v_i > 0\}_{i=1}^{n}$.

Code

Data generation

The data is generated as given in the question.

SeedRandom[1256]

fSurvivalGompertzDistRand[α_, β_] := 
 ProbabilityDistribution[(1/((E^(α/β) Gamma[
          0, α/β])/β) E^(((1 - 
           E^(t β)) α)/β)), {t, 0, ∞}]

data = RandomVariate[fSurvivalGompertzDistRand[0.016, 0.65], {20000}];

σ = 2.5;
dataNoise = 
  data + RandomVariate[NormalDistribution[0, σ], {20000}];

(Re-)start the process

The results of the maximization step are gathered in the lists signalVals and noiseVals.

SeedRandom[5456]
signalVals = {};
noiseVals = {};

Maximization

Select a sample with "good enough" kurtosis. This is not necessary, simple random sampling would do, but it might help getting better results faster.

hk = 1000;
While[! (10 < Abs[hk] < 40),
 dnSample = RandomSample[dataNoise, 40];
 hk = ExKurtosis[dnSample];
 vars = Array[x, Length[dnSample]];
 ]
hk

Solve the maximization problem:

AbsoluteTiming[
 sol = Maximize[
   Join[
    {Abs[ExKurtosis[dnSample - vars]], Abs[ExKurtosis[vars]] < 0.1,
     Abs[Mean[vars]] < 0.05, 
     Abs[σ - StandardDeviation[vars]] < 0.1, 
     Mean[Map[If[Abs[#] < σ, 1, 0] &, vars]] > 0.66 },
    Map[Abs[#] <= 3.1 σ &, vars],
    Map[# > 0 &, dnSample - vars]
    ], vars]
 ]

(* {79.7731, {0.931781, {x[1] -> 0.488303, x[2] -> 0.0693204, 
   x[3] -> -1.58657, x[4] -> -2.73186, x[5] -> -0.792337, 
   x[6] -> -0.301162, x[7] -> 0.0463628, x[8] -> 0.24009, 
   x[9] -> 2.15609, x[10] -> 0.844921, x[11] -> 0.877771, 
   x[12] -> -0.988591, x[13] -> 0.814648, x[14] -> -1.98969, 
   x[15] -> -0.0298853, x[16] -> -0.189145, x[17] -> 0.850365, 
   x[18] -> 0.521628, x[19] -> -1.80022, x[20] -> 0.607911, 
   x[21] -> 0.0872866, x[22] -> 0.68063, x[23] -> -0.0647998, 
   x[24] -> -2.32211, x[25] -> -2.8472, x[26] -> 1.95862, 
   x[27] -> 1.04585, x[28] -> -1.0081, x[29] -> 1.04367, 
   x[30] -> -0.140025, x[31] -> 1.44755, x[32] -> -0.540915, 
   x[33] -> 0.46877, x[34] -> 2.14427, x[35] -> 0.437988, 
   x[36] -> 0.99062, x[37] -> 0.462472, x[38] -> -0.11133, 
   x[39] -> 0.260179, x[40] -> 1.55722}}} *)

While doing the experiments I stopped the maximization process if I thought it takes too much time (more than ~3 minutes).

Accumulate the results

signalVals = Append[signalVals, dnSample - vars /. sol[[2]]];
noiseVals = Append[noiseVals, vars /. sol[[2]]];

opts = {ImageSize -> Medium, PlotRange -> All};
Grid[{{Histogram[Flatten[signalVals], 20, "Probability", opts, 
    PlotLabel -> "Signal"], 
   Histogram[Flatten[noiseVals], 20, "Probability", opts, 
    PlotLabel -> "Noise"]}}]

Reconstruct CDF and PDF

qs = Range[0, 1, 0.1];
xs = Quantile[Flatten[signalVals], qs]

qCDF = Interpolation[Transpose[{xs, qs}], InterpolationOrder -> 1];

Plot[{qCDF[t], 
  Evaluate@CDF[fSurvivalGompertzDistRand[0.016, 0.65], t]}, {t, 
  Min[xs], Max[xs]}, PlotTheme -> "Detailed", 
 PerformanceGoal -> "Speed"]

Plot[{qCDF'[t], 
  Evaluate@PDF[fSurvivalGompertzDistRand[0.016, 0.65], t]}, {t, 
  Min[xs], Max[xs]}, PlotTheme -> "Detailed", 
 PerformanceGoal -> "Speed"]

Monitoring the process

It is helpful to look at goodness of fit measures in order to evaluate the procedure's results.

PearsonChiSquareTest[Flatten[signalVals], 
 fSurvivalGompertzDistRand[0.016, 0.65]]    
(* Out[1087]= 0.061774 *)

PearsonChiSquareTest[Flatten[noiseVals], 
 NormalDistribution[0, σ]]
(* Out[1089]= 0.18782 *)

PearsonChiSquareTest[#, 
   fSurvivalGompertzDistRand[0.016, 0.65]] & /@ signalVals    
(* Out[1090]= {0.301886, 0.238065, 0.142501, 0.80441} *)

PearsonChiSquareTest[#, NormalDistribution[0, σ]] & /@ noiseVals    
(* Out[1091]= {0.46331, 0.608089, 0.970406, 0.338096} *)

Experimental results

Noise with $\sigma = 2.5$

Using noise as provided in the question and making 4 maximization runs, these are the histograms of the obtained distributions:

enter image description here

Here are the reconstructed CDF and PDF:

enter image description here

Noise with $\sigma = 1$

It seems that better results are obtained with smaller standard deviation of the noise. (As expected.) Again using 4 maximization runs. We can see that the CDF is much better approximated.

These are the histograms of the obtained distributions:

enter image description here

These are the reconstructed CDF and PDF:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.