Generate normally distributed random variable using acceptance-rejection method given trail distribution as cauchy

Hello, As part of my work I am trying to solve this task. However I am not able understand how this can be done. Below is my code which I tried but it never worked. Any help on this is really appreciated. I am new to MATLAB.

I am not even able to calculate the value of constant 'c' and the inverse of cdf.

%Define functions and constant c
syms x 
f= (1/sqrt(2*pi)).*exp(-(x.^2)/2);
g=(1/pi).*(1/(1+x.^2));
r = f./g;
r1 = diff(r,x) == 0;

extreme_points = solve(r1,x);
extreme_values = subs(r, x, extreme_points);
[maxX, maxidx] = max(extreme_values);
best_location = extreme_points(maxidx);
best_value = simplify(maxX, 'steps', 50)

c = eval(best_value)

%plot densities f(x), g(x) and function c*g(x)
figure(1)
t=0:0.01:20;
plot(t,f(t),t,g(t),t,c*g(t))

N=10000;
X=zeros(1,N);
for i=1:N
    while 1
        u=rand;        % Simulating from 
        y=-2*log(1-u); % g(x) 
        u=rand;
        if c*g(y)*u<f(y)
            x=y;
            break
        end
    end
    X(i)=x;
end

figure(2)
%Plot normalized histogram together with density. 
hold off
histogram(X,'Normalization','pdf')
hold on
plot(t,f(t),'LineWidth',3)

Answer 1

[Aditya]

It looks like you are trying to implement a rejection sampling method to sample from a normal distribution using a Cauchy distribution as the proposal distribution. Your approach is mostly correct, but there are a few issues and misunderstandings in your MATLAB code that need to be addressed.

syms x
% Define the target density f(x) and proposal density g(x)
f = (1/sqrt(2*pi)) * exp(-x^2 / 2);
g = (1/pi) * (1 / (1 + x^2));

% Calculate the ratio f/g
r = f / g;

% Find the maximum of the ratio f/g to determine the constant c
r1 = diff(r, x) == 0;
extreme_points = solve(r1, x);
extreme_values = subs(r, x, extreme_points);
[maxX, maxidx] = max(double(extreme_values));
best_value = maxX;
c = double(best_value); % Correctly evaluate c

% Plot densities f(x), g(x), and function c*g(x)
figure(1)
t = -10:0.01:10; % Adjusted range for better visualization
plot(t, double(subs(f, x, t)), t, double(subs(g, x, t)), t, c*double(subs(g, x, t)))
legend('f(x)', 'g(x)', 'c*g(x)')
title('Densities f(x), g(x) and c*g(x)')

% Rejection sampling
N = 10000;
X = zeros(1, N);
for i = 1:N
    while true
        u = rand; % Uniform random number
        % Simulate from g(x) using inverse transform sampling
        y = tan(pi * (u - 0.5));
        u = rand; % Another uniform random number for rejection step
        if c * double(subs(g, x, y)) * u < double(subs(f, x, y))
            x = y;
            break
        end
    end
    X(i) = x;
end

% Plot normalized histogram together with density f(x)
figure(2)
histogram(X, 'Normalization', 'pdf')
hold on
plot(t, double(subs(f, x, t)), 'r', 'LineWidth', 2)
title('Sampled Distribution vs. Target Density f(x)')
legend('Sampled Distribution', 'f(x)')

This corrected code should now properly calculate the constant (c), simulate values from the normal distribution using rejection sampling with a Cauchy proposal distribution, and plot the results. Remember, the success of rejection sampling highly depends on how well the proposal distribution (g(x)) covers the target distribution (f(x)), which is controlled by (c).