学校でのT.Aは、最小二乗フィッティングアルゴリズムの例としてこのコードを私に示しました。最小自乗フィッティングのためのCython/numpyと純粋なnumpy
import numpy as np
#return the coefficients (a0,..aN) of the fit y=a0+a1*x+..an*x^n
#with associated sigma dy
#x,y,dy are all np.arrays with dtype= np.float64
def fit_poly(x,y,dy,n):
V = np.asmatrix(np.diag(dy**2))
M = []
for k in range(n+1):
M.append(x**k)
M = np.asmatrix(M).T
theta = (M.T*V.I*M).I*M.T*V.I*np.asmatrix(y).T
cov_t = (M.T*V.I*M).I
return np.asarray(theta.T)[0], np.asarray(cov_t)
私はcythonを使用してコードを最適化しようとしています。私はこのコード
cimport numpy as np
import numpy as np
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef poly_c(np.ndarray[np.float64_t, ndim=1] x ,
np.ndarray[np.float64_t, ndim=1] y np.ndarray[np.float64_t,ndim=1]dy , np.int n):
cdef np.ndarray[np.float64_t, ndim=2] V, M
V=np.asmatrix(np.diag(dy**2),dtype=np.float64)
M=np.asmatrix([x**k for k in range(n+1)],dtype=np.float64).T
return ((M.T*V.I*M).I*M.T*V.I*(np.asmatrix(y).T))[0],(M.T*V.I*M).I
を得たが、実行時には、私は必ず出力はどこが同じになるように、「アサート」を使用した、両方のプログラムのために同じであるように思われます。私は何が間違っている/間違っていますか?
ありがとうございました。うまくいけば、私を助けることができます。
PS:これは[x**k for k in range(n+1)]
を除き(私はこのプロファイリングを呼び出すことができるかどうかを確認しますが、E/Wではない)
import numpy as np
from polyC import poly_c
from time import time
from pancho_fit import fit_poly
#pancho's the T.A,sup pancho
x=np.arange(1,1000)
x=np.asarray(x,dtype=np.float64)
y=3*x+np.random.random(999)
y=np.asarray(y,dtype=np.float64)
dy=np.array([y.std() for i in range(1,1000)],dtype=np.float64)
t0=time()
a,b=poly_c(x,y,dy,4)
#a,b=fit_poly(x,y,dy,4)
print("time={}s".format(time()-t0))
いくつかの中間結果(たとえばM.T * V.I)を保存することをお勧めします。また、速度と数値精度にすぐれた別の行列を直ちに乗算するときに、行列の逆転を回避する方法がよくあります。 – DavidW