import numpy as np import scipy as sc from scipy.optimize import excitingmixing, linearmixing, diagbroyden, root,fsolve import scipy.optimize

Stdevv = 426.77320952046455 Skewv = -0.10601859653075847 Kurtv=
3.1010848699980573 c = [   0.        ,  421.42467836,   -7.35895799,    1.73349647]

def f(p):
     f1=np.power(c[1],2)+6*c[1]*c[3]+2*np.power(c[2],2)+15*np.power(c[3],2)-np.power(Stdevv,2)
     f2=c[2]*(6*np.power(c[1],2)+8*np.power(c[2],2)+72*c[1]*c[3]+270*np.power(c[3],2))-np.power(Stdevv,3)*Skewv
     f3=60*np.power(c[2],4)+3*np.power(c[1],4)+10395*np.power(c[3],4)+60*np.power(c[1],2)*np.power(c[2],2)+4500*np.power(c[2],2)*np.power(c[3],2)
     +630*np.power(c[1],2)*np.power(c[3],2)+936*c[1]*np.power(c[2],2)*c[3]+3780*c[1]*np.power(c[3],3)+60*np.power(c[1],3)*c[3]-np.power(Stdevv,4)*Kurtv
     return(f1,f2,f3) p = c[1:4]

C = scipy.optimize.linearmixing(f,p, iter=1000) print(C)