# -*- coding: utf-8 -*-
"""
Created on Thu Feb  3 18:47:48 2022

@author: Rene
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def running_average(sample):
  n = len(sample)  
  r_average=np.zeros(n)
  sum = 0
  for i in np.arange(0,n):
      sum = sum+sample[i]
      r_average[i] = sum/(i+1)
  return(r_average)

m = 500

alpha=4
beta=2
mean_gamma=alpha/beta
variance_gamma = alpha/(beta**2)
sigma_gamma = variance_gamma**0.5
random_gamma=np.random.gamma(alpha, 1/beta, m)
r_ave_gamma=running_average(random_gamma)


x_points = np.arange(0,1001)/1000
LB_x = 0
UB_x = 4
x_points = (UB_x-LB_x)*x_points+LB_x
norm_100_pdf = norm.pdf(x_points,mean_gamma,sigma_gamma/(100)**0.5)
norm_300_pdf = norm.pdf(x_points,mean_gamma,sigma_gamma/(300)**0.5)

# Plotting a four panel figure

plt.rc('font', family='serif', size='10')
plt.figure(figsize=(10, 5))

Figure = plt.figure(0)
Figure.set_figwidth(9)
Figure.set_figheight(3.5)
Figure.subplots_adjust(hspace=0.5, wspace=0.25)
Figure.tight_layout() 

# Adding a panel to the figure
Panel = Figure.add_subplot(1,2,1)  
LB_x = 0
UB_x = m
off_set_x = 0.025*(UB_x-LB_x)
x_lims = (LB_x-off_set_x,UB_x+off_set_x)
LB_y = 0
UB_y = 10
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = np.arange(0,m)
the_y = random_gamma
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('k')
plt.ylabel('Random Sample')
title_str = r'$Gamma[4,2]-Sample$'
plt.title(title_str,size = '10')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.scatter(the_x,the_y, color='blue',s=1)
plt.axhline(mean_gamma, lw=2, ls = '--',color = 'black',alpha=0.5)

# Adding a panel to the figure
Panel = Figure.add_subplot(1,2,2)  
LB_x = 0
UB_x = m
off_set_x = 0.025*(UB_x-LB_x)
x_lims = (LB_x-off_set_x,UB_x+off_set_x)
LB_y = 1.5
UB_y = 2.5
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = np.arange(0,m)
the_y = r_ave_gamma
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('k')
plt.ylabel('Running Average')
title_str = r'Average of $Gamma[4,2]$ Sample'
plt.title(title_str,size = '10')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.plot(the_x,the_y, color='red', lw=1,alpha=0.75)

plt.plot(100+30*norm_100_pdf,x_points, color='green', lw=1,ls = '--')
plt.plot(300+30*norm_300_pdf,x_points, color='green', lw=1,ls = '--')

plt.axhline(mean_gamma, lw=2, ls = '--',color = 'black',alpha=0.5)
plt.axvline(100, lw=1, ls = ':',color = 'black',alpha=0.5)
plt.axvline(300, lw=1, ls = ':',color = 'black',alpha=0.5)

plt.scatter([100,300],[r_ave_gamma[100],r_ave_gamma[300]],color="blue")


# adding a title to the plot
text_str = 'Graphical Depiction of LOLN & CLT - Average of Gamma[4,2] RV\'s'
plt.suptitle(text_str,size='12',y=0.99)

plt.savefig('Refinement_Law_of_large_numbers.png', dpi=300)