# code example for Metropolis-Hastings MCMC samplig of a Gaussian Distribution
#  by    ChatGPT v. 3.5

# Question to chatGPT: 
# "I am interested in a sequence of numbers sampled from a gaussian
#   distribution using markov chain Monte Carlo. Provide example code.

# this resulted in sub-optimal code with a strong bias towares large values

# 2nd Question:
# I tried the code; the resulting distribution is biased towards high values.
# Can you fix this  in the example ?

# here is the 2nd result (my corrections marked as # !GQ)

# Unfortunately, the problem was not solved (wrong acceptance criterion)
# (see my fix below !).

# a third attempt (pointing to the wrong calculation of the acceptance
# criterion) led to an updated, but still non-correct code !

   # !!! this code by chatGPT is not correct !!!
   # -> chatGPT_GQ_testMCMC.py for corrected version 

import numpy as np
import matplotlib.pyplot as plt

# Define the target Gaussian distribution
target_mean = 5.0
target_stddev = 2.0

# Define the proposal distribution (e.g., another Gaussian)
proposal_mean = 0.0
proposal_stddev = 1.0

# Number of samples to generate
num_samples = 1000

# Initialize the Markov chain
current_sample = np.random.rand()  # Start with an arbitrary initial sample

# Lists to store the generated samples
samples = []

# Metropolis-Hastings algorithm
for _ in range(num_samples):
    # Propose a new sample from the proposal distribution
    proposed_sample = np.random.normal(current_sample, proposal_stddev)

    # Calculate the unnormalized likelihood for the current and proposed samples
    likelihood_current = (
        1 / (np.sqrt(2 * np.pi) * target_stddev)
        * np.exp(-0.5 * ((current_sample - target_mean) / target_stddev) ** 2)
    )
    likelihood_proposed = (
        1 / (np.sqrt(2 * np.pi) * target_stddev)
        * np.exp(-0.5 * ((proposed_sample - target_mean) / target_stddev) ** 2)
    )

    # Calculate the unnormalized proposal likelihood for the current and proposed samples
    proposal_current = (
        1 / (np.sqrt(2 * np.pi) * proposal_stddev)
        * np.exp(-0.5 * ((current_sample - proposal_mean) / proposal_stddev) ** 2)
    )
    proposal_proposed = (
        1 / (np.sqrt(2 * np.pi) * proposal_stddev)
        * np.exp(-0.5 * ((proposed_sample - proposal_mean) / proposal_stddev) ** 2)
    )

    # Calculate the acceptance ratio (Metropolis-Hastings ratio)
    acceptance_ratio = (likelihood_proposed * proposal_current) / (likelihood_current * proposal_proposed)

    # Accept or reject the proposed sample
    if np.random.uniform(0, 1) < min(1, acceptance_ratio):
        current_sample = proposed_sample

    samples.append(current_sample)

# Plot the generated samples
plt.hist(samples, bins=30, density=True, alpha=0.5, color="blue", label="MCMC Samples")
x = np.linspace(-5, 15, 100)
y = (1 / (np.sqrt(2 * np.pi) * target_stddev)) * np.exp(
    -0.5 * ((x - target_mean) / target_stddev) ** 2
)
plt.plot(x, y, color="red", label="True Gaussian")
plt.legend()
plt.xlabel("Sampled Values")
plt.ylabel("Probability Density")
plt.title("Sampling from Gaussian Distribution using MCMC")
plt.show()