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).
