이 답변은 this 게시물과 관련이 있습니다. 여기서 나는 x과 y 오류로 적합하다고 설명합니다. 따라서 ODR 모듈이 필요 없지만 수동으로 수행 할 수 있습니다. 따라서 leastsq 또는 minimize을 사용할 수 있습니다. 제약에 관해서, 나는 가능한 한 그들을 피하려고 노력하는 다른 지위에서 명확히했다. 프로그래밍과 수학에 대한 세부 사항은 약간 성가 시지만, 특히 안정적이고 빠를 것으로 예상되는 경우에는 여기에서 수행 할 수 있습니다. 나는 단지 거친 아이디어를 줄 것이다. 우리가 y0 > m * x0**(-c)을 원한다고 말하십시오. 로그 형식에서는 eta0 > mu - c * xeta0으로 작성할 수 있습니다. 나는. 과 같은 eta0 = mu - c * xeta0 + alpha**2이 있습니다. 다른 불평등에 대해서도 마찬가지입니다. 두 번째 상한은 beta**2이지만 어느 것이 더 작은 지 결정할 수 있으므로 자동으로 다른 조건을 충족시킵니다. 같은 한계가 gamma**2 및 delta**2 인 하한에 적용됩니다. alpha 및 gamma으로 작업 할 수 있습니다. 우리는 불평등 조건을 결합하여이 둘을 관련시킬 수 있습니다. 마지막으로 sigma 및 alpha = sqrt(s-t)* sigma/sqrt(sigma**2 + 1)에 적합 할 수 있습니다. s 및 t은 부등식에서 파생됩니다. sigma/sqrt(sigma**2 + 1) 함수는 alpha이 특정 범위에서 달라 지도록하는 하나의 옵션입니다. 즉, alpha**2 < s-t radicand가 음수가 될 수 있다는 사실은 해가없는 경우를 나타냅니다. alpha을 알고 있으면 mu 및 따라서 m이 계산됩니다. 따라서 적합 매개 변수는 c 및 sigma이며, 이는 부등식을 고려하여 m을 사용합니다. 나는 그것을 지치고 작동하지만 버전은 가장 안정적인 버전이 아닙니다. 요청시 게시 할 수 있습니다.
이미 수작업 잔여 기능이 있으므로 두 번째 옵션이 있습니다. 우리는 단지 우리 자신의 chi**2 함수를 소개하고 제약을 허용하는 minimize을 사용합니다.minimize 및 constraints 키워드 솔루션은 매우 유연하며 잔여 기능은 다른 기능을 위해 쉽게 수정되며 m * x**(-c)은 물론 전체 구조가 매우 유연합니다. 다음과 같습니다 :
import matplotlib.pyplot as plt
import numpy as np
from random import random, seed
from scipy.optimize import minimize,leastsq
seed(7563)
fig1 = plt.figure(1)
###for gaussion distributed errors
def boxmuller(x0,sigma):
u1=random()
u2=random()
ll=np.sqrt(-2*np.log(u1))
z0=ll*np.cos(2*np.pi*u2)
z1=ll*np.cos(2*np.pi*u2)
return sigma*z0+x0, sigma*z1+x0
###for plotting ellipses
def ell_data(a,b,x0=0,y0=0):
tList=np.linspace(0,2*np.pi,150)
k=float(a)/float(b)
rList=[a/np.sqrt((np.cos(t))**2+(k*np.sin(t))**2) for t in tList]
xyList=np.array([[x0+r*np.cos(t),y0+r*np.sin(t)] for t,r in zip(tList,rList)])
return xyList
###function to fit
def f(x,m,c):
y = abs(m) * abs(x)**(-abs(c))
#~ print y,x,m,c
return y
###how to rescale the ellipse to make fitfunction a tangent
def elliptic_rescale(x, m, c, x0, y0, sa, sb):
#~ print "e,r",x,m,c
y=f(x, m, c)
#~ print "e,r",y
r=np.sqrt((x - x0)**2 + (y - y0)**2)
kappa=float(sa)/float(sb)
tau=np.arctan2(y - y0, x - x0)
new_a=r*np.sqrt(np.cos(tau)**2 + (kappa * np.sin(tau))**2)
return new_a
###residual function to calculate chi-square
def residuals(parameters,dataPoint):#data point is (x,y,sx,sy)
m, c = parameters
#~ print "m c", m, c
theData = np.array(dataPoint)
best_t_List=[]
for i in range(len(dataPoint)):
x, y, sx, sy = dataPoint[i][0], dataPoint[i][1], dataPoint[i][2], dataPoint[i][3]
#~ print "x, y, sx, sy",x, y, sx, sy
###getthe point on the graph where it is tangent to an error-ellipse
ed_fit = minimize(elliptic_rescale, x , args = (m, c, x, y, sx, sy))
best_t = ed_fit['x'][0]
best_t_List += [best_t]
#~ exit(0)
best_y_List=[ f(t, m, c) for t in best_t_List ]
##weighted distance not squared yet, as this is done by scipy.optimize.leastsq
wighted_dx_List = [ (x_b - x_f)/sx for x_b, x_f, sx in zip(best_t_List,theData[:,0], theData[:,2]) ]
wighted_dy_List = [ (x_b - x_f)/sx for x_b, x_f, sx in zip(best_y_List,theData[:,1], theData[:,3]) ]
return wighted_dx_List + wighted_dy_List
def chi2(params, pnts):
r = np.array(residuals(params, pnts))
s = sum([ x**2 for x in r] )
#~ print params,s,r
return s
def myUpperIneq(params,pnt):
m, c = params
x,y=pnt
return y - f(x, m, c)
def myLowerIneq(params,pnt):
m, c = params
x,y=pnt
return f(x, m, c) - y
###to create some test data
def test_data(m,c, xList,const_sx,rel_sx,const_sy,rel_sy):
yList=[f(x,m,c) for x in xList]
xErrList=[ boxmuller(x,const_sx+x*rel_sx)[0] for x in xList]
yErrList=[ boxmuller(y,const_sy+y*rel_sy)[0] for y in yList]
return xErrList,yErrList
###some start values
mm_0=2.3511
expo_0=.3588
csx,rsx=.01,.07
csy,rsy=.04,.09,
limitingPoints=dict()
limitingPoints[0]=np.array([[.2,5.4],[.5,5.0],[5.1,.9],[5.7,.9]])
limitingPoints[1]=np.array([[.2,5.4],[.5,5.0],[5.1,1.5],[5.7,1.2]])
limitingPoints[2]=np.array([[.2,3.4],[.5,5.0],[5.1,1.1],[5.7,1.2]])
limitingPoints[3]=np.array([[.2,3.4],[.5,5.0],[5.1,1.7],[5.7,1.2]])
####some data
xThData=np.linspace(.2,5,15)
yThData=[ f(x, mm_0, expo_0) for x in xThData]
#~ ###some noisy data
xNoiseData,yNoiseData=test_data(mm_0, expo_0, xThData, csx,rsx, csy,rsy)
xGuessdError=[csx+rsx*x for x in xNoiseData]
yGuessdError=[csy+rsy*y for y in yNoiseData]
for testing in range(4):
###Now fitting with limits
zipData=zip(xNoiseData,yNoiseData, xGuessdError, yGuessdError)
estimate = [ 2.4, .3 ]
con0={'type': 'ineq', 'fun': myUpperIneq, 'args': (limitingPoints[testing][0],)}
con1={'type': 'ineq', 'fun': myUpperIneq, 'args': (limitingPoints[testing][1],)}
con2={'type': 'ineq', 'fun': myLowerIneq, 'args': (limitingPoints[testing][2],)}
con3={'type': 'ineq', 'fun': myLowerIneq, 'args': (limitingPoints[testing][3],)}
myResult = minimize(chi2 , estimate , args=(zipData,), constraints=[ con0, con1, con2, con3 ] )
print "############"
print myResult
###plot that
ax=fig1.add_subplot(4,2,2*testing+1)
ax.plot(xThData,yThData)
ax.errorbar(xNoiseData,yNoiseData, xerr=xGuessdError, yerr=yGuessdError, fmt='none',ecolor='r')
testX = np.linspace(.2,6,25)
testY = np.fromiter((f(x, myResult.x[0], myResult.x[1]) for x in testX), np.float)
bx=fig1.add_subplot(4,2,2*testing+2)
bx.plot(xThData,yThData)
bx.errorbar(xNoiseData,yNoiseData, xerr=xGuessdError, yerr=yGuessdError, fmt='none',ecolor='r')
ax.plot(limitingPoints[testing][:,0],limitingPoints[testing][:,1],marker='x', linestyle='')
bx.plot(limitingPoints[testing][:,0],limitingPoints[testing][:,1],marker='x', linestyle='')
ax.plot(testX, testY, linestyle='--')
bx.plot(testX, testY, linestyle='--')
bx.set_xscale('log')
bx.set_yscale('log')
plt.show()
제공하는 결과 내가 네 가지 제한 사항 (행) 확인 
############
status: 0
success: True
njev: 8
nfev: 36
fun: 13.782127248002116
x: array([ 2.15043226, 0.35646436])
message: 'Optimization terminated successfully.'
jac: array([-0.00377715, 0.00350225, 0. ])
nit: 8
############
status: 0
success: True
njev: 7
nfev: 32
fun: 41.372277637885716
x: array([ 2.19005695, 0.23229378])
message: 'Optimization terminated successfully.'
jac: array([ 123.95069313, -442.27114677, 0. ])
nit: 7
############
status: 0
success: True
njev: 5
nfev: 23
fun: 15.946621924326545
x: array([ 2.06146362, 0.31089065])
message: 'Optimization terminated successfully.'
jac: array([-14.39131606, -65.44189298, 0. ])
nit: 5
############
status: 0
success: True
njev: 7
nfev: 34
fun: 88.306027468763432
x: array([ 2.16834392, 0.14935514])
message: 'Optimization terminated successfully.'
jac: array([ 224.11848736, -791.75553417, 0. ])
nit: 7
. 결과는 정상적으로 그리고 로그 스케일 (열)로 표시됩니다. 추가 작업을하면 오류가 발생할 수 있습니다.
비대칭 오류에 대한 업데이트 솔직히 말해서이 속성을 처리하는 방법을 모르겠습니다. Naively, 나는 this post과 비슷한 내 비대칭 손실 함수를 정의 할 것입니다. x 및 y 오류가 발생하는 경우, 긍정적이거나 부정적인 측면을 검사하는 대신 4 분면으로 오류를 발생시킵니다. 따라서 내 오류 타원은 연결된 4 개의 조각으로 바뀝니다. 그럼에도 불구하고, 그것은 다소 합리적입니다. 테스트와 작동 방식을 보여주기 위해 선형 함수로 예제를 만들었습니다. OP가 자신의 요구 사항에 따라 두 가지 코드를 결합 할 수 있다고 생각합니다.
선형의 경우는 다음과 같다 적합 :
는
import matplotlib.pyplot as plt
import numpy as np
from random import random, seed
from scipy.optimize import minimize,leastsq
#~ seed(7563)
fig1 = plt.figure(1)
ax=fig1.add_subplot(2,1,1)
bx=fig1.add_subplot(2,1,2)
###function to fit, here only linear for testing.
def f(x,m,y0):
y = m * x +y0
return y
###for gaussion distributed errors
def boxmuller(x0,sigma):
u1=random()
u2=random()
ll=np.sqrt(-2*np.log(u1))
z0=ll*np.cos(2*np.pi*u2)
z1=ll*np.cos(2*np.pi*u2)
return sigma*z0+x0, sigma*z1+x0
###for plotting ellipse quadrants
def ell_data(aN,aP,bN,bP,x0=0,y0=0):
tPPList=np.linspace(0, 0.5 * np.pi, 50)
kPP=float(aP)/float(bP)
rPPList=[aP/np.sqrt((np.cos(t))**2+(kPP*np.sin(t))**2) for t in tPPList]
tNPList=np.linspace(0.5 * np.pi, 1.0 * np.pi, 50)
kNP=float(aN)/float(bP)
rNPList=[aN/np.sqrt((np.cos(t))**2+(kNP*np.sin(t))**2) for t in tNPList]
tNNList=np.linspace(1.0 * np.pi, 1.5 * np.pi, 50)
kNN=float(aN)/float(bN)
rNNList=[aN/np.sqrt((np.cos(t))**2+(kNN*np.sin(t))**2) for t in tNNList]
tPNList = np.linspace(1.5 * np.pi, 2.0 * np.pi, 50)
kPN = float(aP)/float(bN)
rPNList = [aP/np.sqrt((np.cos(t))**2+(kPN*np.sin(t))**2) for t in tPNList]
tList = np.concatenate([ tPPList, tNPList, tNNList, tPNList])
rList = rPPList + rNPList+ rNNList + rPNList
xyList=np.array([[x0+r*np.cos(t),y0+r*np.sin(t)] for t,r in zip(tList,rList)])
return xyList
###how to rescale the ellipse to touch fitfunction at point (x,y)
def elliptic_rescale_asymmetric(x, m, c, x0, y0, saN, saP, sbN, sbP , getQuadrant=False):
y=f(x, m, c)
###distance to function
r=np.sqrt((x - x0)**2 + (y - y0)**2)
###angle to function
tau=np.arctan2(y - y0, x - x0)
quadrant=0
if tau >0:
if tau < 0.5 * np.pi: ## PP
kappa=float(saP)/float(sbP)
quadrant=1
else:
kappa=float(saN)/float(sbP)
quadrant=2
else:
if tau < -0.5 * np.pi: ## PP
kappa=float(saN)/float(sbN)
quadrant=3
else:
kappa=float(saP)/float(sbN)
quadrant=4
new_a=r*np.sqrt(np.cos(tau)**2 + (kappa * np.sin(tau))**2)
if quadrant == 1 or quadrant == 4:
rel_a=new_a/saP
else:
rel_a=new_a/saN
if getQuadrant:
return rel_a, quadrant, tau
else:
return rel_a
### residual function to calculate chi-square
def residuals(parameters,dataPoint):#data point is (x,y,sxN,sxP,syN,syP)
m, c = parameters
theData = np.array(dataPoint)
bestTList=[]
qqList=[]
weightedDistanceList = []
for i in range(len(dataPoint)):
x, y, sxN, sxP, syN, syP = dataPoint[i][0], dataPoint[i][1], dataPoint[i][2], dataPoint[i][3], dataPoint[i][4], dataPoint[i][5]
### get the point on the graph where it is tangent to an error-ellipse
### i.e. smallest ellipse touching the graph
edFit = minimize( elliptic_rescale_asymmetric, x , args = (m, c, x, y, sxN, sxP, syN, syP))
bestT = edFit['x'][0]
bestTList += [ bestT ]
bestA,qq , tau= elliptic_rescale_asymmetric(bestT, m, c , x, y, aN, aP, bN, bP , True)
qqList += [ qq ]
bestYList=[ f(t, m, c) for t in bestTList ]
### weighted distance not squared yet, as this is done by scipy.optimize.leastsq or manual chi2 function
for counter in range(len(dataPoint)):
xb=bestTList[counter]
xf=dataPoint[counter][0]
yb=bestYList[counter]
yf=dataPoint[counter][1]
quadrant=qqList[counter]
if quadrant == 1:
sx, sy = sxP, syP
elif quadrant == 2:
sx, sy = sxN, syP
elif quadrant == 3:
sx, sy = sxN, syN
elif quadrant == 4:
sx, sy = sxP, syN
else:
assert 0
weightedDistanceList += [ (xb - xf)/sx, (yb - yf)/sy ]
return weightedDistanceList
def chi2(params, pnts):
r = np.array(residuals(params, pnts))
s = np.fromiter((x**2 for x in r), np.float).sum()
return s
####...to make data with asymmetric error (fixed); for testing only
def noisy_data(xList,m0,y0, sxN,sxP,syN,syP):
yList=[ f(x, m0, y0) for x in xList]
gNList=[boxmuller(0,1)[0] for dummy in range(len(xList))]
xerrList=[]
for x,err in zip(xList,gNList):
if err < 0:
xerrList += [ sxP * err + x ]
else:
xerrList += [ sxN * err + x ]
gNList=[boxmuller(0,1)[0] for dummy in range(len(xList))]
yerrList=[]
for y,err in zip(yList,gNList):
if err < 0:
yerrList += [ syP * err + y ]
else:
yerrList += [ syN * err + y ]
return xerrList, yerrList
###some start values
m0=1.3511
y0=-2.2
aN, aP, bN, bP=.2,.5, 0.9, 1.6
#### some data
xThData=np.linspace(.2,5,15)
yThData=[ f(x, m0, y0) for x in xThData]
xThData0=np.linspace(-1.2,7,3)
yThData0=[ f(x, m0, y0) for x in xThData0]
### some noisy data
xErrList,yErrList = noisy_data(xThData, m0, y0, aN, aP, bN, bP)
###...and the fit
dataToFit=zip(xErrList,yErrList, len(xThData)*[aN], len(xThData)*[aP], len(xThData)*[bN], len(xThData)*[bP])
fitResult = minimize(chi2, (m0,y0) , args=(dataToFit,))
fittedM, fittedY=fitResult.x
yThDataF=[ f(x, fittedM, fittedY) for x in xThData0]
### plot that
for cx in [ax,bx]:
cx.plot([-2,7], [f(x, m0, y0) for x in [-2,7]])
ax.errorbar(xErrList,yErrList, xerr=[ len(xThData)*[aN],len(xThData)*[aP] ], yerr=[ len(xThData)*[bN],len(xThData)*[bP] ], fmt='ro')
for x,y in zip(xErrList,yErrList)[:]:
xEllList,yEllList = zip(*ell_data(aN,aP,bN,bP,x,y))
ax.plot(xEllList,yEllList ,c='#808080')
### rescaled
### ...as well as a scaled version that touches the original graph. This gives the error shortest distance to that graph
ed_fit = minimize(elliptic_rescale_asymmetric, 0 ,args=(m0, y0, x, y, aN, aP, bN, bP))
best_t = ed_fit['x'][0]
best_a,qq , tau= elliptic_rescale_asymmetric(best_t, m0, y0 , x, y, aN, aP, bN, bP , True)
xEllList,yEllList = zip(*ell_data(aN * best_a, aP * best_a, bN * best_a, bP * best_a, x, y))
ax.plot(xEllList, yEllList, c='#4040a0')
###plot the fit
bx.plot(xThData0,yThDataF)
bx.errorbar(xErrList,yErrList, xerr=[ len(xThData)*[aN],len(xThData)*[aP] ], yerr=[ len(xThData)*[bN],len(xThData)*[bP] ], fmt='ro')
for x,y in zip(xErrList,yErrList)[:]:
xEllList,yEllList = zip(*ell_data(aN,aP,bN,bP,x,y))
bx.plot(xEllList,yEllList ,c='#808080')
####rescaled
####...as well as a scaled version that touches the original graph. This gives the error shortest distance to that graph
ed_fit = minimize(elliptic_rescale_asymmetric, 0 ,args=(fittedM, fittedY, x, y, aN, aP, bN, bP))
best_t = ed_fit['x'][0]
#~ print best_t
best_a,qq , tau= elliptic_rescale_asymmetric(best_t, fittedM, fittedY , x, y, aN, aP, bN, bP , True)
xEllList,yEllList = zip(*ell_data(aN * best_a, aP * best_a, bN * best_a, bP * best_a, x, y))
bx.plot(xEllList, yEllList, c='#4040a0')
plt.show()
상부 그래프는 기존의 선형 함수를 나타내고, 일부 데이터가 비대칭 가우시안 에러를 이용하여 이로부터 생성
나타내는한다. 오류 막대가 그려지며, 조각 별 오류 타원 (회색 ... 및 선형 함수를 만지려면 재조정 됨, 파란색). 아래쪽 그래프는 fit 함수뿐만 아니라 재조정 된 piecewise ellipses를 보여줍니다.
오 ... 단지 주목하고 있습니다, 상한선입니까, 아니면 하한선도 있습니까? ... 여하튼, 내 해결책은 양쪽 모두를 허용합니다. – mikuszefski
@mikuszefski 예, 상한선입니다 만, y와 x 에러가 비대칭이 될 필요가 있음을 언급하는 것을 잊었습니다. 즉, 아래쪽과 위쪽의 에러에 대해 yerr_l과 yerr_u를 넣어야 할 필요가 있습니다. 엑스! – astrostudent
오, 알았어. 그 패키지를 제대로 가지고 있다면, ODR로도 작동하지 않을 것이다. 실제로, 직교 거리가 비대칭 오류에 대한 것이 아닌지 확실하지 않습니다. 아마도 사분면적일 수 있습니다. 하지만 'x'와 'y'는 상관 관계가 없습니까? – mikuszefski