import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# GAMMA

mean = 6000  # Mean (mu)
sd = 3000  # Standard deviation (sigma)

# Calculate the shape (k) and scale (theta) parameters
k = (mean ** 2) / (sd ** 2)
theta = (sd ** 2) / mean

# Draw 1000 samples from the gamma distribution
samples = np.random.gamma(k, theta, 1000)

# Calculate the mean and standard deviation of the samples
sample_mean = np.mean(samples)
sample_std = np.std(samples)

print(sample_mean)
print(sample_std)

clipped_samples = np.clip(samples, 4000, 15000)
clipped_sample_mean = np.mean(clipped_samples)
clipped_sample_std = np.std(clipped_samples)
print(clipped_sample_mean)
print(clipped_sample_std)

#plt.hist(clipped_samples, bins=50, color='skyblue', edgecolor='black', alpha=0.7)
sns.histplot(data=clipped_samples, kde=True)
plt.title('Histogram of 1000 Samples from Gamma Distribution')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plt.show()
"""

# HALFNORMAL
scale = 5000
# Draw 1000 samples from a normal distribution and take their absolute values
samples_halfnormal = np.abs(np.random.normal(0, scale, 4000))
mean_halfnormal = np.mean(samples_halfnormal)
std_halfnormal = np.std(samples_halfnormal)
print(mean_halfnormal)
print(std_halfnormal)

clipped_samples = np.clip(samples_halfnormal, 4000, 15000)
clipped_sample_mean = np.mean(clipped_samples)
clipped_sample_std = np.std(clipped_samples)

plt.hist(clipped_samples, bins=50, color='salmon', edgecolor='black', alpha=0.7)
plt.title('Histogram of 1000 Samples from Half-Normal Distribution')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plt.show()
"""