import numpy as np
import matplotlib.pyplot as plt
filename = 'E:/data.txt'

x = np.loadtxt(filename)[:,0]
y = np.loadtxt(filename)[:,1]
#plt.plot(x,y,'*')
b0 = [9.38,-0.03,25348.9123,0.05,4.0]
for i in range(5):
    f0 = b0[0]
    A = b0[1]
    x0 = b0[2]
    sigma = b0[3]
    gamma = b0[4]
    g = 1.0-np.exp(-np.power(x-x0,2)/2./sigma**2) 
    f = f0+A*np.power(g,gamma)   
    #plt.plot(x,y,'b*')
    #plt.plot(x,f,'r*')
    #plt.show()
    X = np.ones((len(x),5))
    X[:,1] = np.power(g,gamma)    #derivacia podla A
    X[:,2] = A*gamma*np.multiply(np.multiply(np.power(g,gamma-1),g-1),x-x0)/sigma**2   #derivacia podla x0
    X[:,3] = A*gamma*np.multiply(np.multiply(np.power(g,gamma-1),g-1),np.power(x-x0,2))/sigma**3   #derivacia podla sigma
    X[:,4] = A*np.power(g,gamma)*np.log(g)    #derivacia podla gamma
    dY=y-f
    V = np.dot(np.transpose(X),X)
    U = np.dot(np.transpose(X),dY)
    b0 += np.dot(np.linalg.inv(V),U)
    s =  np.sqrt(np.dot(dY.T,dY)/(len(X)-len(X[0])))
    db1 = s*np.sqrt(np.diag(np.linalg.inv(V)))

b = []
for j in range(1000):
    x = np.loadtxt(filename)[:,0]
    y = np.loadtxt(filename)[:,1]
    ind = np.random.choice(range(len(x)),len(x))
    x = x[ind]
    y = y[ind]
    b0 = [9.38,-0.03,25348.9123,0.05,4.0]
    for i in range(5):
        f0 = b0[0]
        A = b0[1]
        x0 = b0[2]
        sigma = b0[3]
        gamma = b0[4]
        g = 1.0-np.exp(-np.power(x-x0,2)/2./sigma**2) 
        f = f0+A*np.power(g,gamma)   
        X = np.ones((len(x),5))
        X[:,1] = np.power(g,gamma)    #derivacia podla A
        X[:,2] = A*gamma*np.multiply(np.multiply(np.power(g,gamma-1),g-1),x-x0)/sigma**2   #derivacia podla x0
        X[:,3] = A*gamma*np.multiply(np.multiply(np.power(g,gamma-1),g-1),np.power(x-x0,2))/sigma**3   #derivacia podla sigma
        X[:,4] = A*np.power(g,gamma)*np.log(g)    #derivacia podla gamma
        dY=y-f
        V = np.dot(np.transpose(X),X)
        U = np.dot(np.transpose(X),dY)
        b0 += np.dot(np.linalg.inv(V),U)
        s = np.sqrt(np.dot(dY.T,dY)/(len(X)-len(X[0])))
    b.append(b0)

b = np.array(b)
print(db1)
print(np.std(b,axis = 0))