Convert symbolic to numeric code: speed up morphing

Numeric solution

In this example I compute the $W_2$ geodesic (Wasserstein distance) between two densities defined as InterpolatingFunction.

(* unnormalized density functions *)
uf = Interpolation[{{-2, .5}, {0, 2}, {.5, 1}, {1, .5}}];
ug = Interpolation[{{-1, 1}, {0, .5}, {1, 2}, {2, .5}}];

(* normalized density functions *)
f[x_] = uf[x]/NIntegrate[uf[x], {x, -2, 1}];
g[x_] = ug[x]/NIntegrate[ug[x], {x, -1, 2}];
ℱ = ProbabilityDistribution[f[x], {x, -2, 1}];
\[ScriptCapitalG] = ProbabilityDistribution[g[x], {x, -1, 2}];

Show[
 Plot[f[x], {x, -2, 1}, PlotStyle -> Blue, Filling -> 0],
 Plot[g[x], {x, -1, 2}, PlotStyle -> Red, Filling -> 0],
 PlotRange -> {All, {0, All}}, AxesOrigin -> {0, 0}]

The points xF are a linear sampling of the domain of f. The points qF are the quantiles associated to the points xG. The points xℱ are the union of the two, in order to ensure that both densities are discretized sufficiently well.

xF = Range[-2, 1, .05];
xG = Range[-1, 2, .05];
qF = InverseCDF[ℱ, CDF[\[ScriptCapitalG], xG]];
qG = InverseCDF[\[ScriptCapitalG], CDF[ℱ, xF]];
xℱ = Union[xF, qF];
x\[ScriptCapitalG] = Union[xG, qG];

X[t] is the interpolation between the starting and final points, whereas dens[t] is the intermediate density at those points.

X[t_] := (1 - t) xℱ + t x\[ScriptCapitalG]
dens[t_] := 1/((1 - t)/f /@ xℱ + t/g /@ x\[ScriptCapitalG])

The resulting density can be visualized as

ListLinePlot[Evaluate@Table[{X[t], dens[t]}\[Transpose], {t, 0, 1, .1}]]

The transport map can also be computed and plotted with

dT = f /@ xℱ/g /@ x\[ScriptCapitalG];
T = Interpolation[{{xℱ}\[Transpose], x\[ScriptCapitalG], dT}\[Transpose]];
Plot[T[x], {x, xℱ[[1]], xℱ[[-1]]}]

Symbolic solution

Mathematica appears to be able to deal with distributions, CDF, inverse CDF and pushforwards of distributions:

    ℱ = UniformDistribution[-1 + {-1, 1}/2];
    \[ScriptCapitalG] = TriangularDistribution[1 + {-1, 1}];
    T[x_] = InverseCDF[\[ScriptCapitalG], CDF[ℱ, x]] // Simplify;
    \[ScriptCapitalD][t_] := TransformedDistribution[(1 - t) x + t T[x], x \[Distributed] ℱ]
    Plot[{PDF[ℱ, x], PDF[\[ScriptCapitalG], x]}, {x, -2, 2}]
    Plot[Evaluate@Table[PDF[\[ScriptCapitalD][t], x], {t, 0., 1., .1}], {x, -2, 3}]

Symbolic integration of UnitBox and UnitTriangle

While it's true that

    Integrate[UnitBox[y], {y, -∞, x}]

and

    Integrate[UnitTriangle[y], {y, -∞, x}]

do not work as intended, giving a slight hint regarding the domain of x helps in both cases

    Integrate[UnitBox[y], {y, -∞, x}, Assumptions -> x ∈ Reals]
    Integrate[UnitTriangle[y], {y, -∞, x}, Assumptions -> x ∈ Reals]

and the returned result are piecewise functions. An antiderivative can also be found with

    Derivative[-1][UnitBox][x]
    Derivative[-1][UnitTriangle][x]

You can use NDSolveValue to create an interpolating function representation of the inverse. Basically, suppose you want to invert f. Then:

f[finv[x]] == x

where finv is the inverse function. So, an ODE for the inverse function is:

D[f[finv[x]] == x, x]

f'[finv[x]] finv'[x] == 1

Let's use this for your G function:

g[x_] := UnitTriangle[x-3]
G[x_] := Integrate[g[s], {s, -Infinity, x}]

Then we have:

Ginv = NDSolveValue[{G'[inv[x]] inv'[x] == 1, inv[G[3]] == 3}, inv, {x, 0, 1}]

However, it's easy to see that we can use g instead of G', so it will be quicker to do:

Ginv = Quiet @ NDSolveValue[{g[inv[x]] inv'[x] == 1, inv[G[3]] == 3}, inv, {x, 0, 1}];

The quieted messages are associated with the the fact that g is zero when x is at one of the endpoints, 0 or 1. Let's check:

G[Ginv[0]]
G[Ginv[.5]]
G[Ginv[.75]]
G[Ginv[1]]

0.

0.5

0.75

1.

So, Ginv is an interpolating function representation of the inverse of G, and you can take derivatives of it as desired, e.g.:

D[Ginv[Sin[x]], x] /. x->3

-1.86349