# Plot a transformed distribution with an If function

Suppose I have a normal distribution $A \sim \mathcal{N}(\mu,\sigma^{2})$ with a known cutoff point (percentile) on this distribution called $c$. Based on this point the distribution needs to be transformed which leads to a new distribution called $W$. Next, suppose $s$ equals a deterministic value. The first part of $A \in [-\infty,c]$ becomes equal to a transformation of $A$ more specifically $\mathcal{N}(\mu+s,\sigma^{2})$, the second part $]c,+\infty]$ remains equal to $A$ To be clear, $c$ and $s$ are known in advance and do not need to be estimated.

To summarise: $W = \begin{cases} A+s & \text{if } A \leq c \\ A & \text{if } A > c \end{cases}$

I used the following coded to obtain the mean and variance of $W$:

a = Exp[-(x - μ)^2/(2 σ^2)]/(Sqrt[2 Pi] σ);
domain[a] = {x, -Infinity, Infinity} && {Element[μ, Reals], σ > 0, Element[c, Reals], s > 0};
w = If[x <= c, x + s, x];
meand = Expect[w, a]
vard = Var[w, a]


Is there a quick way to plot the pdf of $W$ using Mathematica ? (I also have mathStatica)

• Have you tried defining your function with Piecewise and assigning some values to $s$, $\mu$, $\sigma$? – Aisamu Nov 18 '14 at 10:41

This is actually a very nice question, and although it appears deceptively simple, it is not.

The problem

Let $$X \sim N(\mu, \sigma^2)$$ with pdf $$f(x)$$, and let $$s>0$$ and $$c$$ denote known constants. That is:

Find the pdf of $$W$$, where:

$$W = \begin{cases} X+s & \text{if } X \leq c \\ X & \text{if } X > c \end{cases}$$

The problem is not as simple as it looks, because the transformation from $$X$$ to $$W$$ is NOT one to one. In particular, there is a zone of overlap where both branches of the piecewise function can contribute to the same value of $$W$$.

Solution

The cdf of $$W$$ is simply $$P(W:

where:

1. the Prob function is from the mathStatica package for Mathematica, which the OP is using;
2. the True zone denotes $$c.

The pdf of $$W$$ is the derivative of the cdf wrt $$w$$:

with domain of support:

 domain[pdf] = {w, -Infinity, Infinity}  && {Element[c, Reals], s > 0};


All done.

Plotting the pdf of $$W$$

Given parameters, say:

params = {μ -> 2, σ -> 1, s -> 3, c -> 2};


... here is a plot of the pdf of $$W$$:

PlotDensity[pdf /. params, {w, -1, 7}]


Monte Carlo check

Here are some different parameters ($$c$$ has changed):

params = {μ -> 2, σ -> 1, s -> 3, c -> 2.8};


Here are $$500,000$$ pseudo-random drawings from $$X$$ and $$W$$:

xdata = RandomReal[NormalDistribution[2, 1], {500000}];
wdata = Map[If[# < 2.8, # + 3, #] &, xdata];


Now we can compare:

• the theoretical pdf (dashed RED curve) derived above
• to the empirical pdf (squiggly BLUE curve)

...

 FrequencyPlot[wdata, {0, 7, .02}, pdf /. params]


Looks fine. :)

If you don't have mathStatica ...

If you don't have mathStatica, one can still obtain the transformed pdf using just Mathematica with:

 mmasol = PDF[TransformedDistribution[Piecewise[{{x + s, x <= c}, {x, x > c}}],
Distributed[x, NormalDistribution[μ, σ]],
Assumptions -> {Element[c, Reals], s > 0}], w]


... but the solution produced is not very natural, and is expressed in terms of unnecessary DiracDelta functions, UnitStep functions and Erf functions:

(1/2)*DiracDelta[c + s - w]*(1 + Erf[(w - μ)/(Sqrt[2]*σ)] -
(Erf[(c - μ)/(Sqrt[2]*σ)] - Erf[(w - μ)/(Sqrt[2]*σ)] -
Erfc[(s - w + μ)/(Sqrt[2]*σ)])*(-1 +  UnitStep[c - w])) +(1/2)*DiracDelta[   c - w]*(-Erfc[(s - w + μ)/(Sqrt[2]*σ)] +   (-Erf[(c - μ)/(Sqrt[2]*σ)] +    Erf[(w - μ)/(Sqrt[2]*σ)] +
Erfc[(s - w + μ)/(Sqrt[2]*σ)])*(1 -
UnitStep[-c - s + w])) + (1/(Sqrt[2*Pi]*σ))*(UnitStep[c - w]/
E^((s - w + μ)^2/(2*σ^2)) +
(E^(-((w - μ)^2/(2*σ^2))) +
E^(-((s - w + μ)^2/(2*σ^2))))*(1 -
UnitStep[c - w])*    (1 - UnitStep[-c - s + w]) +
UnitStep[-c - s + w]/E^((w - μ)^2/(2*σ^2)))


# Important Update!

This answer is wrong, for reasons that, even after @wolfies explanation, I had a hard time grasping! I'll leave it here as warning just as an example of not what to do! (It is mathematically unsound, to say the least)

A = Exp[-(x - μ)^2/(2 σ^2)]/(Sqrt[2 Pi] σ);
A1 = A /. (μ -> μ + s);
W = Piecewise[{{A1, x <= c}}, A];

Manipulate[
Plot[{W} /. {μ -> mu, σ -> sig, s -> ss, c -> cut}, {x, -10, 10}, PlotRange -> {{-10, 10}, {0, .5}}]
, {{mu, 1}, -10, 10}, {{sig, 2}, -10, 10}, {{ss, 2}, -10, 10}, {{cut, 1}, -10, 10}]


• How do you draw the vertical line joining the two pieces of Piecewise? Usually I just get a discontinuous plot. I'm studying your code but I don't see it. – becko Nov 27 '14 at 16:37
• That was a mistake, since there was supposed to be an exclusion there. You can force that behavior with ExclusionsStyle -> Automatic! "PerformanceGoal"->"Speed" also does it, but for the wrong reasons... – Aisamu Nov 27 '14 at 17:02

Using the same approach as in @wolfies' answer (i.e., differentiating the CDF instead of using the PDF), gives the same result without the DiracDeltas:

tdist = TransformedDistribution[If[u <= c, u + s, u],
u \[Distributed] NormalDistribution[μ, σ], Assumptions -> {Element[{c, s}, Reals]}]
(* or Boole[u <= c] s + u instead of If[u <= c, u + s, u] *)

FullSimplify[CDF[tdist, w]]


pdF[ww_, cc_, ss_, μ0_, σ0_] := FullSimplify[D[CDF[tdist, w], w] /.
{μ -> μ0, σ -> σ0,  s -> ss, c -> cc}];

pdF[w, c, s, μ, σ]


If you add the assumption s>0, you get the same result as in @wolfies' answer:

pdF[ww_, cc_, ss_, μ0_, σ0_] := FullSimplify[D[CDF[tdist, w], w] /.
{μ -> μ0, σ -> σ0, s -> ss, c -> cc}, Assumptions -> {Element[c, Reals], s > 0}];
pdF[w, c, s, μ, σ]


Plot[Evaluate[pdF[w, 2, 3, 2, 1.]], {w, -1, 7}, BaseStyle -> Thick, Filling -> Axis]