import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm,uniform,t

nx=1000
ny=1000

X = norm(0.58,1)
Y = uniform() # Y ~ U[0,1]

x = X.rvs(nx)
y = Y.rvs(ny)

# plt.figure()
# plt.plot(x,'o')
# plt.savefig('x_sample.png')

# plt.figure()
# plt.plot(y,'o')
# plt.savefig('y_sample.png')

bar_x = x.mean()
bar_y = y.mean()

var_x = np.sum((x - bar_x)**2)/(nx-1)
var_y = np.sum((y - bar_y)**2)/(ny-1)


var_bar_x = var_x/nx
var_bar_y = var_y/ny

t_obs = (bar_x-bar_y)/(np.sqrt(var_bar_x + var_bar_y))

print('t_obs=',t_obs)

alpha = 0.05

p_value = 1 - t.cdf(np.abs(t_obs),df=nx+ny-2) + t.cdf(-np.abs(t_obs),df=nx+ny-2)

print('p_value=', p_value)

if p_value < alpha:
    print('Reject Ho')
else:
    print('Do not reject Ho')