# -*- coding: utf-8 -*-
"""
Created on Fri Oct 11 00:04:32 2024

@author: jotzu
"""

import numpy as np
import matplotlib.pyplot as plt
sigma = 2
NList = [2,3,4,5,7,10,15,20,35,50]
mMeanList = []
vMeanList = []
vUnbiasedMeanList = []
sMeanList = []
sUnbiasedMeanList = []
sUnbiased2MeanList = []
for N in NList:    
    M=100000
    a = np.random.normal(0,sigma,size=(N,M))
    m = np.mean(a,axis=0)
    v = np.var(a,axis=0)
    s = np.sqrt(v)
    vUnbiased = np.var(a,axis=0)*(N/(N-1))
    sUnbiased = np.sqrt(vUnbiased)
    sUnbiased2 = np.std(a,axis=0)*np.sqrt((N/(N-1.5)))
    mMeanList.append(m.mean())
    vMeanList.append(v.mean())
    vUnbiasedMeanList.append(vUnbiased.mean())
    sMeanList.append(s.mean())
    sUnbiasedMeanList.append(sUnbiased.mean())
    sUnbiased2MeanList.append(sUnbiased2.mean())
print(a)
plt.plot(NList,mMeanList,marker="o")
plt.hlines(0,0,max(NList),color="grey")
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated mean")
plt.show()
plt.close()
#%%

plt.plot(NList,vMeanList,marker="o")
plt.hlines(sigma**2,0,max(NList),color="grey")
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated variance")
plt.show()
plt.close()


#%%

plt.plot(NList,vMeanList,marker="o")
plt.hlines(sigma**2,0,max(NList),color="grey")
n = np.arange(2,max(NList),.01)
plt.plot(n,4*(n-1)/(n),color="red")
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated variance")
plt.show()
plt.close()


#%%

plt.plot(NList,vMeanList,marker="o")
plt.plot(NList,vUnbiasedMeanList,marker="o")
plt.hlines(sigma**2,0,max(NList),color="grey")
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated variance")
plt.show()
plt.close()
#%%

plt.plot(NList,vUnbiasedMeanList,marker="o",color="orange")
plt.hlines(sigma**2,0,max(NList),color="grey")
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated variance")
plt.show()
plt.close()


#%%

plt.plot(NList,sMeanList,marker="o")
plt.plot(NList,sUnbiasedMeanList,marker="o")
plt.hlines(sigma,0,max(NList),color="grey")
# n = np.arange(2,max(NList),.01)
# plt.plot(n,1*n/(n-1))
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated std. dev.")

plt.show()
plt.close()
#%%

plt.plot(NList,sMeanList,marker="o")
plt.plot(NList,sUnbiasedMeanList,marker="o")
plt.plot(NList,sUnbiased2MeanList,marker="o")
plt.hlines(sigma,0,max(NList),color="grey")
# n = np.arange(2,max(NList),.01)
# plt.plot(n,1*n/(n-1))
plt.xlabel("Size of Dataset")
plt.ylabel("Estimated std. dev.")
plt.show()
plt.close()