Как решить нелинейную функцию для вектора параметров в python - PullRequest
Новая
       73

Как решить нелинейную функцию для вектора параметров в python

0 голосов
/ 20 апреля 2020

У меня есть нелинейная функция, которую я хочу решить при различных значениях параметров.

Вот MWE: 

import numpy  as np 
import tensorflow as tf
import scipy.optimize 

# fixed parameters
kon = 0.01
mu  = 1.5
fi  = 0.5 
kappa = 22
w = 0.63

# varying parameters 
n =100
xs = tf.random.normal(shape=(n,), stddev=0.2)
eps = tf.random.normal(shape=(n,), stddev=0.17)
z = tf.sigmoid(tf.random.normal(shape=(n,), stddev=0.22))

Мое решение 

def get_leisure(z, eps, x0):
    hvec = np.empty((0,))
    # leisure today
    for ze,ei,xs in zip(z, eps, x0):
        ei=np.exp(ei)
        xs=np.exp(xs)
        # define the function for leisure 
        def leisure_function(hi):
            return (mu/fi)*np.log(hi) -(1-mu)*kappa*(hi)**(1+(1/fi))- mu*(np.log(w*ei*xs)-np.log(kon))-np.log(ze)

        htemp = scipy.optimize.newton_krylov(leisure_function, 0.5)
        hvec = np.append(hvec, htemp)
    return hvec

Мой вопрос:

Поскольку число дел, которое мне нужно l oop, чтобы решить для hi неизвестно, может быть большим, есть ли способ сделать это лучше? например, чтобы избежать l oop?

Я не опытный пользователь python, и я был бы признателен за любую помощь.

1 Ответ

1 голос
/ 21 апреля 2020

Просто чтобы проиллюстрировать то, что я хотел бы отметить в своем комментарии.

Если я запустил ваш код таким, какой он есть, и рассчитал время: 

import numpy  as np 
import tensorflow as tf
import scipy.optimize 

from timeit import default_timer as timer
# fixed parameters
kon = 0.01
mu  = 1.5
fi  = 0.5 
kappa = 22
w = 0.63

# varying parameters 
n=10000
x0 = tf.random.normal(shape=(n,), stddev=0.2)
eps = tf.random.normal(shape=(n,), stddev=0.17)
z = tf.sigmoid(tf.random.normal(shape=(n,), stddev=0.22))

def get_leisure(z, eps, x0):
    hvec = np.empty((0,))
    # leisure today
    for ze,ei,xs in zip(z, eps, x0):
        ei=np.exp(ei)
        xs=np.exp(xs)
        # define the function for leisure 
        def leisure_function(hi):
            return (mu/fi)*np.log(hi) -(1-mu)*kappa*(hi)**(1+(1/fi))- mu*(np.log(w*ei*xs)-np.log(kon))-np.log(ze)

        htemp = scipy.optimize.newton_krylov(leisure_function, 0.5)
        hvec = np.append(hvec, htemp)
    return hvec

start = timer()
msh855_result = get_leisure(z, eps, x0)
end = timer()
print(f'elapsed time: { end - start} s')

Это займет около 30.04823809700065 s в универсальный c станок.

Та же проблема в той же машине, но с векторизованным подходом 

start = timer()
e_eps = np.exp(eps)
e_x0 = np.exp(x0)
c_1=mu/fi * np.ones(n)
c_2=(1-mu)*kappa * np.ones(n)
c_3= 1 + (1/fi)
c_4=mu*np.log(w*e_eps*e_x0/kon)+np.log(z)
def fun(x):
    return c_1[0] * np.log(x) - c_2[0] * ((x) ** c_3)-c_4
v_result = scipy.optimize.newton_krylov(fun, 0.5 * np.ones(n))
end = timer()
print(f'elapsed time: { end - start} s')

Это займет всего 0.051494838000508025 s

И если вы ' Вы не уверены, что результат довольно близок: 

>>> sum(np.sqrt(((msh855_result - v_result)**2)/n))
2.0087031341897715e-06
...