Extracting signal from gaussian noise

Question

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:

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:

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.

Answer 1

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:

Here are the reconstructed CDF and PDF:

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:

These are the reconstructed CDF and PDF: