from scipy.optimize import curve_fit
import numpy as np

def func_model(x,a,b):
    return a*x[0]+b*x[1]


def calc_hara(val,f,pi):
    xdata=np.array([val,f],dtype="float")
    pi = np.array(pi,dtype="float")
    # print(xdata.shape)
    estimates, vcov=curve_fit(func_model,xdata,pi)
    # print(estimates)
    gamma=1-(.235-.03)/(.734**2)/estimates[0]
    eta=estimates[1]*(.734)**2/(.235-.03)

    return gamma,eta

gammas=[]
etas=[]
for j in range(1000):
    val=[1]
    for k in range(10):
        val.append(np.random.normal(val[-1],.05))
    val=val[1:]
    pi= np.random.randint(5,95,size=10)/100
    f=[]
    for i in range(5):
        f.append(np.exp(-.03*(4-i)))
    f.extend(f)

    gamma,eta = calc_hara(val,f,pi)
    gammas.append(gamma)
    etas.append(eta)

    print(gamma, eta, val[-1])

print(min(gammas),max(gammas),np.mean(gammas))
print(min(etas),max(etas),np.mean(etas))