2
$\begingroup$

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)

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

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))) 
| improve this answer | |
$\endgroup$
1
$\begingroup$

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

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – becko Nov 27 '14 at 16:37
  • $\begingroup$ 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... $\endgroup$ – Aisamu Nov 27 '14 at 17:02
1
$\begingroup$

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

enter image description here

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

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

enter image description here

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, μ, σ]

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.