#!/usr/bin/env nix-shell
#!nix-shell -i python3 -p python3 python3Packages.pandas python3Packages.numpy python3Packages.scipy python3Packages.pyarrow
# coding: utf-8

from scipy.integrate import quad
import pandas as pd
import numpy as np

from pathlib import Path


def normalProbabilityDensity(x):
    constant = 1.0 / np.sqrt(2*np.pi)
    return (constant * np.exp((-x**2) / 2.0))

home = Path.home()
snt_path = Path(home, ".standard-normal.csv")

if snt_path.exists():
    snt = pd.read_csv(snt_path, index_col=0)
else:
    snt = pd.DataFrame(data    = [],
                       index   = np.round(np.arange(0, 3.5, 0.1), 2),
                       columns = np.round(np.arange(0.00, 0.1, 0.01), 2))

    for index in snt.index:
        for column in snt.columns:
            z = np.round(index + column, 2)
            value, _ = quad(normalProbabilityDensity, np.NINF, z)
            value = np.round(value - 0.5, 4)
            snt.loc[index,column] = value
    snt.to_csv(snt_path)

def find_z(p):
    posns = []
    res = snt.isin([p])
    srs = res.any()
    cols = list(srs[srs == True].index)
    for col in cols:
        rows = list(res[col][res[col] == True].index)
        for row in rows:
            posns.append(float(col) + row)
    return np.max(posns)

mean = float(input("Sample Mean: "))
std = float(input("Standard Deviation: "))
N = int(input("Sample Size: "))
alpha = float(input("Alpha [0, 1): "))
alpha_2 = alpha / 2

Z_alpha_2 = find_z(alpha_2)

plus_minus = np.round(Z_alpha_2 * (std / np.sqrt(N)), 4)

print(f"{mean} +/- {plus_minus}")
print(f"[{mean - plus_minus}, {mean + plus_minus}]")

### Local Variables:
### mode: python
### End: