struggling to solve problem

I have this problem to solve:

Implementing Customized Probability Distributions [50 pts]

Frequently in engineering, you have to draw random samples from a customized probability distribu-

tion function p(x). MATLAB provides functions such as rand and randn that draw samples from

uniform and Normal distributions, respectively. But, what if your probability function is something

else?

(a) In this problem, you will implement a method to draw samples from:

p(x) =

(− 2

25 x + 2

5 if x ∈ [0, 5],

0 otherwise. (1)

To do this, first explicitly calculate the cumulative probability density function:

P (x) =

Z x

−∞

p(u)du. (2)

Note: You do not have to use MATLAB for this part. Show your result in your report.

(b) The range of all cumulative probability density functions is [0, 1] because of the requirement

that R ∞

−∞ p(x)dx = 1. With this fact, you can sample p(x) by first drawing a random sample

y from a uniform distribution on the domain [0, 1] (this is what rand does) and then get x by

using x = P −1(y). Create a function called myRand that provides a random sample x from the

probability density function p(x).

function x = myRand()

% Outputs:

% x: the random sample drawn from p(x)

(c) Draw 100,000 samples from your myRand function and create a histogram of the results us-

ing MATLAB ’s histogram function. Set the Normalization property in histogram to

pdf. Refer to https://www.mathworks.com/help/matlab/ref/matlab.graphics.

chart.primitive.histogram.html to look up how to do this. Also, on the histogram,

plot p(x) in a different color. Screenshots are not accepted. Comment on how well your im-

plementation approximates the probability density function., e.g., ”Fig. 1 shows Y quantity as a

function of X quantity. We observe that [a few sentences on the findings].”

I solved part a then wrote the following code to solve the rest:

%% Part b
function x = myRand()
y = rand();  % Generate a random number between 0 and 1
if y <= 0  % Handle the case when y is 0
    x = 0;
elseif y >= 1  % Handle the case when y is 1
    x = 5;
else
    x = (-5/2) + sqrt(25*y + 25/4);  % Invert the CDF equation
end
%% Part c
% Generate 100,000 samples
samples = zeros(1, 100000);
for i = 1:100000
    samples(i) = myRand();
end

% Create the histogram
figure;
histogram(samples, 'Normalization', 'pdf');
hold on;

% Plot the original probability density function p(x)
x = linspace(0, 5, 1000);
p = (-2/25)*x + 2/5;
plot(x, p, 'r', 'LineWidth', 2);

% Add labels and title
xlabel('x');
ylabel('Probability Density');
title('Histogram of Random Samples and p(x)');

% Show the plot
grid on;
legend('Samples', 'p(x)');
end

HOWEVER, instead of getting a histogram that follows the overall trend of p(x) plotted too, i got a histogram that basically goes the opposite way (inserted below) Can yoiu help me. thank you.

[Walter Roberson]

It may interest you to know that none of the floating point random number generation techniques supported by MATLAB can generate exactly 0 or exactly 1 (unless you want to count generation on the normal distribution.) So when you use rand() you know that the value will be strictly > 0 and < 1 -- your else is the only clause that can get invoked.

This also has the implication that you can vectorize the random number generation.

Answer 1

[Torsten]

Your inverted CDF is wrong.

In myRand, set x to

x = 5*(1-sqrt(1-y))

instead of

x = (-5/2) + sqrt(25*y + 25/4);

And - as Walter said - you don't need this if-statement because y will always satisfy 0 < y < 1.