# -*- coding: utf-8 -*-
"""
Created on Tue Aug 31 18:37:42 2021

@author: Rene
"""

import numpy as np
import matplotlib.pyplot as plt

n = 50
p = np.arange(1,n)/n
mu = -np.log(p)

MSE_S_Theo  =p*(1-p)/30
MSE_S_values = np.zeros(n-1)
MSE_T_values = np.zeros(n-1)

# Generate m random samples from Poisson distribution and evaluate S and T Estimates
m = 10000
i_range = np.arange(0,m)
j_range = np.arange(0,n-1)

for j in j_range:
    S = np.zeros(m)
    T = np.zeros(m)
    for i in i_range:
        x = np.random.poisson(mu[j],30)
        mu_hat=np.mean(x)
        y_hat=(30-np.count_nonzero(x))/30
        S[i]=y_hat
        T[i]=np.exp(-mu_hat)
    MSE_S_values[j]=np.mean((S-p[j])**2)
    MSE_T_values[j]=np.mean((T-p[j])**2)   

# Plotting the analysis results    
LB = 0
UB = 1  
off_set_x = 0.02
off_set_y = 0.00025
max_y = 0.01
    
# Setting the attributes for the Figure     
plt.rc('font', family='serif', size='10')
plt.figure(figsize=(10, 5))
MSE_Figure = plt.figure()
MSE_Figure.subplots_adjust(hspace=0.4, wspace=0.4)
x_lims = (LB-off_set_x,UB+off_set_x)
y_lims = (0-off_set_y,max_y+off_set_y)
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)

x_values = [0.1,0.25,0.5,0.75,0.9]
for xc in x_values:
    plt.axvline(x=xc,ls=':',lw=0.5,color='black')

plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('p')
plt.ylabel('MSE')

plt.plot(p,MSE_S_Theo,ls='--',lw=0.5,label='MSE S Theo',color='black')
plt.plot(p,MSE_S_values,ls='-',lw=0.5,label='MSE S Emp',color='red')
plt.plot(p,MSE_T_values,ls='-',lw=0.5,label='MSE T Emp',color='blue')

plt.legend(loc = 'best')
MSE_Figure.suptitle('MSE for Estimators S and T',size='14')
plt.savefig('S_T_MSE_p.png', dpi=1200)