# Disentangling the data

An experiment yields a functional dependence on two variables $F(x,y)$, can be imagined as a 2D map. It is known that the function can be factored as follows $F(x,y)=fx(x) gy(y) +gx(x) fy(y)$. Given the tabulated $F(x,y)$ I would like to numerically find the four unknown functions $fx$, $fy$, $gx$, $gy$. It is clear that in full generality the problem is ill posed and has no unique solution. There are however several constraints that are known a priori:

i) $fx(x)$, $fy(y)$ are non-zero only on a finite interval

ii) $gy(y)$ has a gaussian profile, $gx(x)$ decay exponentially for $x>0$.

Just to give an example let us use the following definitions:

fx[x_] := 1/(Exp[-x - 7] + 1) + 1/(Exp[x - 7] + 1) - 1
gy[y_] := Exp[-y^2]
fy[y_] := (1/(Exp[-y - 10] + 1) + 1/(Exp[y - 10] + 1) - 1) (5 + 0.5 Sin[y])
gx[x_] := Exp[-0.4 x]/(Exp[-5 x] + 5)
F[x_, y_] := fx[x] gy[y] + gx[x] fy[y]


looking as The resulting function has a cross shape

 Plot3D[F[x, y], {x, -15, 15}, {y, -15, 15}, PlotRange -> All, PlotPoints -> {50, 50}] I do not have many ideas how to find $fx,fy,gx,gy$ numerically. I was trying to use a fitting procedure, however, it works like you get what you put in. I am seeking for a more general solution, perhaps with the help of Fourier or wavelet analysis. I would be grateful for any ideas.

• there is a typo here gx[x_] := gx[x_] May 1 '14 at 18:10
• you could use BSplines to parametrize the rest? as in mathematica.stackexchange.com/a/10997/1089 and then use maximum likelihood on the data (which I assume is made of pairs of {x_i,y_i}s ? May 1 '14 at 18:12
• @b.gatessucks and @ chris: thank you for correcting a typo May 2 '14 at 7:43
• Just to make the problem statement more precise: experimental data looks like the 3d plot above, i.e. a set of points $x_,y_i, F(x_i,y_i)$; based on this input the 4 single-argument functions (like on the Fig. 1) need to be numerically reconstructed. May 2 '14 at 7:46
• A general solution might be quite difficult. Therefore I fell a need to formulate 2 simpler problems: i) numerically find coordinates of the intersection point; ii) numerically verify if the angle between 2 waves is strictly 90 degrees. May 2 '14 at 7:53

So, you have a function $F(x,y) = f_x(x)g_y(y) + g_x(x)f_y(y)$, and you want to recover $f_x,g_y,g_x,f_y$.

If you've tabulated the values of $F(x,y)$ in a matrix $\mathbf F$ with entries $f_{ij} = F(x_i,y_j)$, then this amounts to decomposing the matrix as $$\mathbf F \approx \mathbf f_x\mathbf g_y^T + \mathbf g_x\mathbf f_y^T,$$ where $\mathbf f_x,\mathbf g_y,\mathbf g_x,\mathbf f_y$ are column vectors with non-negative entries. (I'm using $\approx$ instead of $=$ because your data presumably has some noise in it.) Stick the vectors into two matrices, $\mathbf W = \begin{bmatrix}\mathbf f_x & \mathbf g_x\end{bmatrix}$ and $\mathbf H = \begin{bmatrix}\mathbf g_y^T \\ \mathbf f_y^T\end{bmatrix}$, and you have $$\underbrace{\mathbf F}_{n\times n} \approx \underbrace{\mathbf W}_{n\times2}\underbrace{\mathbf H}_{2\times n}$$ where all three matrices have non-negative entries. This is precisely the problem of non-negative matrix factorization. And look, there's a Mathematica implementation in the open source Mathematica for Prediction project.

Let's try it!

f = Table[F[x, y], {x, -15, 15, 0.1}, {y, -15, 15, 0.1}];
Needs["NonNegativeMatrixFactorization"];
{w, h} = GDCLS[f, 2];
fx = w[[All, 1]];
gx = w[[All, 2]];
gy = h[[1, All]] // Normal;
fy = h[[2, All]] // Normal;
ListLinePlot[{fx, gx}, PlotRange -> All]
ListLinePlot[{gy, fy}, PlotRange -> All]  There's a bit of cross-talk between the components, and the results appear to be a little different every time you run it (maybe because of random initialization), but overall it looks pretty good.

I am not sure where your problem lies exactly. If you have a set of points, x_i,y_i which obey the PDF F[x,y], you could do maximum likelihood analysis. A parametric model could be

gxa[x_, a_] = Exp[-a x]/(Exp[-5 x] + 5);
G[x_, y_, a_] = fx[x] gy[y] + gxa[x, a] fy[y];


with the corresponding normalization (so that its a PDF)

norma[a_] =
Table[{a,
NIntegrate[G[x, y, a], {x, -20, 0, 20}, {y, -20, 0, 20},
PrecisionGoal -> 2]}, {a, 0.1, 0.9, 0.1}] //
Interpolation[#, a] &;


So that the various PDFs look like this

Table[ContourPlot[G[x, y, a]/norma[a], {x, -8, 8}, {y, -8, 8},
PlotRange -> All, PlotPoints -> {50, 50},
MeshFunctions -> Function[{x, y, z}, z],
ContourStyle -> ColorData[a*10], ContourShading -> False],
{a, 0.1, 0.9, 0.2}] // Show A draw would be drawn from and set of points

sample = RandomVariate[UniformDistribution[{-15, 15}], {10000, 3}];
sample = sample /. {x_, y_, z_} :> {x, y, Abs[z]/10/norma[0.4]};


using the 'keep below PDF' prescription (for a=0.4)

ok = Select[sample, #[] < G[#[], #[], 0.4]/norma[0.4] &];
ok = ok /. {x_, y_?NumberQ, _} -> {x, y}; ok // Length


so that we can check that the maximum likelihood value of our draw corresponds to a=0.4

Table[{a, Plus @@ Log@Map[(G[#[], #[], a]/norma[a]) &, ok]}, {a,
0.1, 0.9, 0.025}] // ListLinePlot which indeed peaks near 0.4.

Note that the above parametric model could be made more general, using e.g. BSplines at the expense of a more complex optimization problem.

Note finally that mathematica has its own MaximumLikelihood function.

• thank you for your answer. It provides a partial solution and demonstrates how difficult a general solution might be. Your solution shows how to find a single parameter when 3 other functions are known as well as the functional dependence on the parameter $a$. I still hope that all 4 functions can be reconstructed by imposing some constraints. May 2 '14 at 7:51

Keep in mind this is only a partial solution and it is supposed to give you ideas and some insight.

First to give you an idea of where I am going consider the following example

img = Import["http://i.stack.imgur.com/yV8FW.png"]
WaveletImagePlot[DiscreteWaveletTransform[img]] Now on to the more interesting problem at hand

fx[x_] := 1/(Exp[-x - 7] + 1) + 1/(Exp[x - 7] + 1) - 1
gy[y_] := Exp[-y^2]
fy[y_] := (1/(Exp[-y - 10] + 1) + 1/(Exp[y - 10] + 1) - 1) (5 + 0.5 Sin[y])
gx[x_] := Exp[-0.4 x]/(Exp[-5 x] + 5)
F[x_, y_] := fx[x] gy[y] + gx[x] fy[y]

dwd = DiscreteWaveletTransform[Table[F[x, y], {x, -15, 15}, {y, -15, 15}]]

WaveletMatrixPlot[dwd] Now do you see why I started with an example ?

Now we shall recover the waves

ListPlot3D[InverseWaveletTransform[
WaveletMapIndexed[#1 0.0 &, dwd, {___, 1 | 4}]], PlotRange -> All] ListPlot3D[InverseWaveletTransform[
WaveletMapIndexed[#1 0.0 &, dwd, {___, 2 | 4}]], PlotRange -> All] Or directly from the wavelet transform

ListPlot3D[Abs@Reverse@dwd[{2}, "Values"], PlotRange -> All,
Boxed -> False, ImageSize -> 500] ListPlot3D[Abs@Reverse@dwd[{1}, "Values"], PlotRange -> All,
Boxed -> False, ImageSize -> 500]
` You can and shall always play with the different parameter settings, use different wavelet families, etc.

• It seems there are many solutions to my problem! May 12 '14 at 15:00