PDF of $Q$ Random Variable

Question

Let $X\sim N(0,25)$, $Y\sim N(10,100)$, $Z\sim N(-10,50)$ and $Q=\tan^{-1}\left(\frac{Z}{\sqrt {X^2+Y^2}}\right)$

When I simulate $Q$ random variable with Monte Carlo method, I'm getting this probability density function graph that you can examine in below.

$Q$ Random Variable Probability Density Function Graph:

And this is also matlab code for the graph.

clc
clear all
close all
for i=1:200000
   x=normrnd(0,25);
   y=normrnd(10,100);
   z=normrnd(-10,50);

   q(i)=atan(z/(sqrt(x*x+y*y)));
end
hist(q,100)
title('Q')

Can I use this graph to guess probability density function of $Q$ random variable? Or what can I do to guess probability density function of $Q$ random variable by using Matlab?

And lastly, what can I do to verify the pdf result that I found in Matlab? How can I determine the reliability of result?

Answer 1

A few thoughts, first when using MATLAB, it would probably pay to use the vectorization capabilities as this would be much faster than the for loop:

N = 1e7;
X = normrnd(0,25,N,1);
Y = normrnd(10,100,N,1);
Z = normrnd(-10,50,N,1);
Q = atan(Z./sqrt(X.^2 + Y.^2));
histogram(Q,'normalization','pdf','binmethod','scott','edgealpha',0)

The option pdf gives you a PDF estimate based on the number of observations in the bin and you can see more details by checking out histogram in MATLAB's documentation.

In terms of figuring out how good it is, I can only recommend increasing $N$. I imagine if you used $N = 10$ or $N = 100$ and compared to the distribution of $N = 1e7$ it would look different. You can use goodness-of-fit test such as Kolmogorov-Smirnoff to compare with different $N$'s if you really desired.