Hej mam do zminimalizowania dość skomplikowaną funkcję która zależy od wielu parametrów i która zawiera sporo funkcji matematycznych typu sinh, cos, pierwiastek itp, czyli funkcje które mają konkretną dziedzine i w momencie gdy chcemy obliczyc jej wartosc w niedozwolonym argumencie to dostajemy NaN. To samo dotyczy dzielenia przez 0 oraz gdy chcemy obliczyc bardzo duza potege jakiejs liczby. Używam funkcji brute z pakietu scipy.optimize i w tej funkcji podaje zakres wartosc dla kazdego argumentu oraz krok o jaki kazdy parametr ma sie zmieniac. W rezultacie dla niektorych wartosci argumentow dostaje wartosc minimlizowanej funkcji NaN i caly algorytm sie psuje i zamiast zwrocic najmniejsza wartosc sposrod wszystkich sprawdzanych to dostaje NaN. Co zrobic zeby NaN nie mialy wplywu na dzialanie algorytmu?
Wydaje mi się, że sporo by pomogło gdybyś jednak pokazał tą funkcję.
katakrowa napisał(a):
Wydaje mi się, że sporo by pomogło gdybyś jednak pokazał tą funkcję.
Ok oto kod:
def cf(u, T, r, alpha, beta, delta):
''' Funkcja charakterystyczna w modelu Kou '''
b = r - q - 2 * delta * (np.log(np.cos(beta/2))-np.log(np.cos((alpha+beta)/2)))
value = np.exp(1j * u * b * T)*(np.cos(beta/2)/np.cosh((alpha*u-1j*beta)/2))**(2*delta*T)
return value
N=2048
L=10
# This class defines puts and calls
class OptionType(enum.Enum):
CALL = 1.0
PUT = -1.0
def CallPutOptionPriceCOSMthd(cf,CP,S0,r,T,K,alpha, beta, delta,N,L):
if K is not np.array:
K = np.array(K).reshape([len(K),1])
# Assigning i=sqrt(-1)
i = np.complex(0.0,1.0)
x0 = np.log(S0 / K)
# Truncation domain
# a = 0.0 - L * np.sqrt(T)
# b = 0.0 + L * np.sqrt(T)
c1=T*alpha*delta*np.tan(beta/2)
c2=0.5*T*alpha**2*delta/((np.cos(beta/2))**2)
# c4=(3*delta*T+2-np.cos(beta))*alpha**4*delta*T/(4*(np.cos(beta/2))**4)
# c4=-1/8*T*alpha**4*delta*(1/8*np.cos(3*beta)-9/8*np.cos(beta)-1)/(np.cos(beta/2))**8
c4=0.25*T*alpha**4*delta*(2-np.cos(beta))/((np.cos(beta/2))**4)
a=c1-L*np.sqrt(c2+np.sqrt(c4))
b=c1+L*np.sqrt(c2+np.sqrt(c4))
# Summation from k = 0 to k=N-1
k = np.linspace(0,N-1,N).reshape([N,1])
u = k * np.pi / (b - a);
# Determine coefficients for put prices
H_k = CallPutCoefficients(CP,a,b,k)
mat = np.exp(i * np.outer((x0 - a) , u))
temp = cf(u, T, r, alpha, beta, delta) * H_k
temp[0] = 0.5 * temp[0]
value = np.exp(-r * T) * K * np.real(mat.dot(temp))
return np.exp(-q*T)*S0+ value - K*np.exp(-r*T)
# Determine coefficients for put prices
def CallPutCoefficients(CP,a,b,k):
if CP==OptionType.CALL:
c = 0.0
d = b
coef = Chi_Psi(a,b,c,d,k)
Chi_k = coef["chi"]
Psi_k = coef["psi"]
if a < b and b < 0.0:
H_k = np.zeros([len(k),1])
else:
H_k = 2.0 / (b - a) * (Chi_k - Psi_k)
elif CP==OptionType.PUT:
c = a
d = 0.0
coef = Chi_Psi(a,b,c,d,k)
Chi_k = coef["chi"]
Psi_k = coef["psi"]
H_k = 2.0 / (b - a) * (- Chi_k + Psi_k)
return H_k
def Chi_Psi(a,b,c,d,k):
psi = np.sin(k * np.pi * (d - a) / (b - a)) - np.sin(k * np.pi * (c - a)/(b - a))
psi[1:] = psi[1:] * (b - a) / (k[1:] * np.pi)
psi[0] = d - c
chi = 1.0 / (1.0 + np.power((k * np.pi / (b - a)) , 2.0))
expr1 = np.cos(k * np.pi * (d - a)/(b - a)) * np.exp(d) - np.cos(k * np.pi
* (c - a) / (b - a)) * np.exp(c)
expr2 = k * np.pi / (b - a) * np.sin(k * np.pi *
(d - a) / (b - a)) - k * np.pi / (b - a) * np.sin(k
* np.pi * (c - a) / (b - a)) * np.exp(c)
chi = chi * (expr1 + expr2)
value = {"chi":chi,"psi":psi }
return value
def Meixner_error_function(p0):
global m, local_opt, opt1
alpha, beta, delta=p0
se = []
for row, option in calls.iterrows():
model_value = CallPutOptionPriceCOSMthd(cf,OptionType.PUT,S0,r,option['T'],(option['strike'],),alpha, beta, delta,N,L)[0][0]
se.append((model_value - option['price']) ** 2)
MSE = sum(se) / len(se)
if m % 1 == 0:
print('%4d |' % m, np.array(p0), '| %7.3f' % MSE)
m += 1
# if local_opt:
# penalty = np.sqrt(np.sum((np.array(p0)-np.array(opt1)) ** 2))
# return MSE + penalty
return MSE**0.5
opt1 = spo.brute(Meixner_error_function, ((0.001,2.51, 0.2),(-np.pi,np.pi, 0.5), (0.01,2,0.16)), finish=None)
Minimalizuje funkcje Meixner_error_function. W nawiasach sa przedzialy dla kazdego argumentu. Niestety funkcja jest na tyle skomplikowana ze po prostu w pewnym momencie liczona jest jej wartosc w nieprawidlowych punktach i w rezultacie dostaje NaN. Najlepszym rozwiazaniem byloby zeby w takim przypadku Python zamiast zwracac NaN zwracał np. 1000000, wowczas jest to na tyle wysoka wartosc ze nie przeszkadzalaby w prawidlowym dzialaniu algorytmu
Funkcje nie mogą Ci zwracać NaN, zdebuguj, żeby zobaczyć co tam się naprawdę dzieje.
No ale dla niektorych argumentow np liczymy e^2000 a taka liczba jest za duża. To samo dotyczy liczenia logarytmu z liczby ujemnej - nie da sie. Chce tak zrobic zeby takie przypadki nie psuly algorytmu a recznie nie da rady ich wykluczyc. Na zdjeciu widac kolejne parametry i wartosc minimalizowanej funkcji dla nich. W pewnym momencie parametry sa takie ze wyskakuje blad i dalej algorytm juz nie bierze pod uwage tych normalnychliczb tylko zwraca parametry dlaktorych funkcja ma wartosc NaN
Debugowałeś? Jak program ma liczyć poprawnie z logarytmem z liczby ujemnej?
Chce tak zrobic zeby takie przypadki nie psuly algorytmu
Obawiam się że tak to się nie da zrobić. Przecież masz tu jakiś algorytm optymalizacyjny. Cała idea takich algorytmów polega na tym, że maja jakiś "kierunek", czy to będzie jakiś gradient descent czy coś innego to bez różnicy. Jak nagle zmienisz wynik na coś losowego, to cała ta optymalizacja idzie w las przecież.
Zresztą "poprawnym" rozwiązaniem byłoby jednak ustalenie dziedziny twojej funkcji dokładnie.
Ad samej precyzji obliczeń, może warto pomyśleć o użyciu gmpy2 i mpfr?