# -*- 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 gamma

the_lambda = 0.5
n = [1,4,9,16]
#n = [25,40,50,100]
the_colors=['blue','red','green','orange']
x = np.arange(0,1001)/1000
LB_x = 0
UB_x = 12
x = (UB_x-LB_x)*x+LB_x
pdf_gamma_1 = gamma.pdf(x,a = n[0],scale = 1/(n[0]*the_lambda))
pdf_gamma_2 = gamma.pdf(x,a = n[1],scale = 1/(n[1]*the_lambda))
pdf_gamma_3 = gamma.pdf(x,a = n[2],scale = 1/(n[2]*the_lambda))
pdf_gamma_4 = gamma.pdf(x,a = n[3],scale = 1/(n[3]*the_lambda))

# 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(4.75)
Figure.subplots_adjust(hspace=0.5, wspace=0.25)

# Adding a panel to the figure
Panel = Figure.add_subplot(2,2,1)  
LB_x = 0
UB_x = 4
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 = 2
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = x
the_y = pdf_gamma_1
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('x')
plt.ylabel('PDF $f(x)$', size = 8)
title_str = r'PDF of average of '+ f'{n[0]:2.0f}'+' RV\'s'+' : $X \sim Exp(0.5)$'
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=the_colors[0], lw=1)
x_values = [2,2]
y_values = [0,gamma.pdf(2,a = n[0],scale = 1/(n[0]*the_lambda))]
plt.plot(x_values,y_values, color='black', lw=1, ls=':')

text_str = r'n = '+ f'{n[0]:2.0f}'
plt.text(0.1,1.8,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
text_str = r'mean = '+ f'{1/the_lambda:2.1f}'
plt.text(0.1,1.6,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
st_dev = (((1/the_lambda)**2)/n[0])**0.5
text_str = r'st.dev. = '+ f'{st_dev:2.1f}'
plt.text(0.1,1.4,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')


# Adding a panel to the figure
Panel = Figure.add_subplot(2,2,2)  
LB_x = 0
UB_x = 4
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 = 2
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = x
the_y = pdf_gamma_2
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('x')
plt.ylabel('PDF $f(x)$', size = 8)
title_str = r'PDF of average of '+ f'{n[1]:2.0f}'+' RV\'s'+' : $X \sim Exp(0.5)$'
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=the_colors[1], lw=1)
x_values = [2,2]
y_values = [0,gamma.pdf(2,a = n[1],scale = 1/(n[1]*the_lambda))]
plt.plot(x_values,y_values, color='black', lw=1, ls=':')

text_str = r'n = '+ f'{n[0]:2.0f}'
plt.text(0.1,1.8,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
text_str = r'mean = '+ f'{1/the_lambda:2.1f}'
plt.text(0.1,1.6,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
st_dev = (((1/the_lambda)**2)/n[1])**0.5
text_str = r'st.dev. = '+ f'{st_dev:2.2f}'
plt.text(0.1,1.4,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')

# Adding a panel to the figure
Panel = Figure.add_subplot(2,2,3)  
LB_x = 0
UB_x = 4
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 = 2
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = x
the_y = pdf_gamma_3
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('x')
plt.ylabel('PDF $f(x)$', size = 8)
title_str = r'PDF of average of '+ f'{n[2]:2.0f}'+' RV\'s'+' : $X \sim Exp(0.5)$'
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=the_colors[2], lw=1)
x_values = [2,2]
y_values = [0,gamma.pdf(2,a = n[2],scale = 1/(n[2]*the_lambda))]
plt.plot(x_values,y_values, color='black', lw=1, ls=':')

text_str = r'n = '+ f'{n[0]:2.0f}'
plt.text(0.1,1.8,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
text_str = r'mean = '+ f'{the_lambda:2.1f}'
plt.text(0.1,1.6,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
st_dev = (((1/the_lambda)**2)/n[2])**0.5
text_str = r'st.dev. = '+ f'{st_dev:2.2f}'
plt.text(0.1,1.4,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')


# Adding a panel to the figure
Panel = Figure.add_subplot(2,2,4)  
LB_x = 0
UB_x = 4
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 = 2
off_set_y = 0.05*(UB_y-LB_y)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

the_x = x
the_y = pdf_gamma_4
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('x')
plt.ylabel('PDF $f(x)$', size = 8)
title_str = r'PDF of average of '+ f'{n[3]:2.0f}'+' RV\'s'+' : $X \sim Exp(0.5)$'
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=the_colors[3], lw=1)
x_values = [2,2]
y_values = [0,gamma.pdf(2,a = n[3],scale = 1/(n[3]*the_lambda))]
plt.plot(x_values,y_values, color='black', lw=1, ls=':')

text_str = r'n = '+ f'{n[3]:2.0f}'
plt.text(0.1,1.8,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
text_str = r'mean = '+ f'{1/the_lambda:2.1f}'
plt.text(0.1,1.6,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')
st_dev = (((1/the_lambda)**2)/n[3])**0.5
text_str = r'st.dev. = '+ f'{st_dev:2.2f}'
plt.text(0.1,1.4,text_str,color = 'red', size = 8,horizontalalignment  = 'left', verticalalignment  = 'bottom')


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