Plot a transformed distribution with an If function

Question

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)

Answer 1

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]

Answer 2

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<w)$:

where:

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

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))) 

Answer 3

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}]