Plotting a transformed distribution as a PDF

Question

I'm rather new to Mathematica having been a longtime Matlab user, so apologies if my question is framed incorrectly or I have missed out expected infomation. I have tried to distill the problem to the essentials.

Effectively I am trying to plot a a transformed distribution in a similar fashion to this example code:

    data = RandomVariate[NormalDistribution[6, 5/3], 10^6];
    Show[Histogram[data, Automatic, "PDF"], 
    Plot[PDF[NormalDistribution[6, 5/3], x], {x, 0, 14}, 
    PlotStyle -> Thick]]

I have create a transformed distribution, sampled some data from it using RandomVariate and plotted it in a histogram.

    hiod = 6.5/2;
    transDistrib = 
    TransformedDistribution[hiod/Tan[theta*Degree], 
    theta \[Distributed] NormalDistribution[6, 5/3]];
    dataTrans = RandomVariate[transDistrib, 10^6];
    Histogram[dataTrans, {0, 80, 1}, "PDF", ImageSize -> Large ]

This works as expected (I have done the same thing in Matlab). However, the equivalent generation of a PDF fails (code just hangs and aborts).

    Plot[PDF[transDistrib, x], {x, 0, 80}, PlotStyle -> Thick]

I get the distinct feeling I an making a terribly simple error, either in my basic coding or understanding of how to represent probability distributions in Mathematica.

Ultimately I would like to derive the equation for the transformed distribution, but I have stumbled at this earlier stage using Mathematicas inbuilt functions.

Any pointers to solving my error would be greatly appreciated. I'm using Mathematica 11.0.1 on macOSX Sierra.

Answer 1

hiod = 6.5/2;

Since the tails of a normal distribution go outside of the allowed range, use TruncatedDistribution to constrain the range of theta

transDistrib = 
  TransformedDistribution[hiod/Tan[theta*Degree], 
   theta \[Distributed] 
    TruncatedDistribution[{0, 90}, NormalDistribution[6, 5/3]]];

dataTrans = RandomVariate[transDistrib, 10^6];

Show[
 Histogram[dataTrans, {0, 80, 1}, "PDF", ImageSize -> Large],
 Plot[PDF[transDistrib, x], {x, 0, 80}, PlotStyle -> Thick]]

EDIT:

With a different mean for the underlying NormalDistribution, the bin specification for the Histogram must be adjusted.

transDistrib = 
  TransformedDistribution[hiod/Tan[theta*Degree], 
   theta \[Distributed] 
    TruncatedDistribution[{0, 90}, NormalDistribution[3, 5/3]]];

dataTrans = RandomVariate[transDistrib, 10^6];

Show[Histogram[dataTrans, {0, 250, 1}, "PDF", ImageSize -> Large], 
 Plot[PDF[transDistrib, x], {x, 0, 250}, PlotStyle -> Thick]]

Answer 2

If most (meaning greater than 99.9%) of the probability is between angles of -90 and +90 degrees, then using the standard transformation of a random variable techniques (find inverse function, multiply by Jacobian, etc.) will give an "exact approximate" or "approximately exact" result. (The exact result requires dealing with a wrapped normal distribution.)

There's no mathematical problem using negative angles and if negative angles are not a problem with your experimental design, the following should work. The density function of $y=hiod/\tan{(\theta \pi/180)}$ where $\theta \sim N(\mu,\sigma)$ (and $\theta$ is in degrees) can be calculated:

μ = 6;
σ = 5/3;
hiod = 13/4;
pdf[y_, hiod_, μ_, σ_] :=
  (2 π σ^2)^(-1/2) Exp[-( (180 ArcCot[y/hiod])/π - μ)^2/(2 σ^2)]*
  180/(hiod π (1 + (y/hiod)^2))

We can take some random samples:

d = NormalDistribution[\[Mu], \[Sigma]];
z = RandomVariate[d, 10^6];
dataTrans = hiod/Tan[z*Degree];

There are many extreme values that make using the defaults for Histogram unsatisfactory so we create a custom set of bins:

lower = 0;
upper = 85;
bins = {Flatten[{Min[dataTrans], Table[i, {i, lower, upper}], Max[dataTrans]}]};
Show[Histogram[dataTrans, bins, "PDF", ImageSize -> Large,
  PlotRange -> {{lower, upper}, Automatic},
  PlotLabel -> Style["μ = 6", Bold, Large]],
 Plot[pdf[y, hiod, μ, σ], {y, lower, upper}]]

Below are the figures for μ=3 and μ=0 showing the histogram and probability density function: