Почему эта реализация оптимизации Nelder Mead не сходится - PullRequest
Новая
       13

Почему эта реализация оптимизации Nelder Mead не сходится

0 голосов
/ 04 октября 2019

Я пытаюсь оценить вероятность успешного схода к оценке максимального правдоподобия для модели логистической регрессии, используя процедуру оптимизации Nelder Mead. Я тестирую различные значения параметров (а именно количество ковариат, а также уровень сигнала), генерирую несколько наборов данных с использованием этих значений параметров и пытаюсь запустить оптимизацию для каждого из них.

Моя проблема заключается в том, чтоНелдер Мид постоянно говорит, что не может сходиться, потому что «максимальное число итераций превышено»;однако выходные параметры кажутся приличными.

Итак, я начинаю задумываться, не требует ли реализация алгоритма необоснованно большого количества итераций, чтобы он сходился / не сходился. У меня ограниченные формальные знания по оптимизации, хотя я подозреваю, что проблема заключается в «размере скачков», которые алгоритм использует, чтобы найти новый симплекс для тестирования и градиент в пространстве параметров. Я прочитал ответ в этой пересекающейся публикации , но, похоже, подразумевается, что нет единого критерия для успешного завершения. Таким образом, я все еще в растерянности из-за того, как вылечить мою проблему в scipy.

Мой код ниже. В данный момент он выводит много значений, поскольку я пытался выяснить, в чем проблема. 

#Script for probability of success of finding maximum likelihood estimator

#Load modules
from scipy.optimize import minimize, root
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import math
import scipy
from scipy.special import erf
from scipy import stats
from matplotlib import cm
import pandas as pd
import matplotlib as mpl
norm = mpl.colors.Normalize(vmin=-1.,vmax=1.)
import time

epsilon = np.finfo(float).eps
df=500


#Define logistic function
def h(beta,X):
    beta = np.reshape(beta, (len(beta),1))
    try:
        a = X @ beta
    except:
        print('beta, X dimensions not aligned')

    h = 1/(1+np.exp(-a))
    return h

def data(p,n,beta):
    #Perhaps put all of the distributions here
    X = scipy.stats.t.rvs(df, size=(n,p))

    #X = scipy.random.normal(0,1, size=(n,p))
    #X = scipy.random.uniform(-5,5, size=(n,p))
    f=h(beta,X)
    y=np.zeros(n)
    for i,val in enumerate(f):
        r=np.random.uniform(0,1)
        if r>val:
            y[i]=1
    return X,y





def successprob(beta0,gamma0,Ncov,n,N,VN):

    Ncov,n,N = np.int(Ncov),np.int(n),np.int(N)
    s=0
    betatries={}
    times=[]
    for sample in range(N):
        beta = np.random.normal(1, 1, Ncov)
        beta = beta * gamma0/VN
        X,y = data(Ncov,n,beta)

        def negloglikelihood(beta_est):
            hlist = h(beta_est,X)
            hlist = np.reshape(hlist, (n,))
            denom = 1-hlist
            hlist = list(hlist)
            denom = list(denom)
            hlist = [val if val>2*epsilon else 2*epsilon for val in hlist]
            denom = [val if val>2*epsilon else 2*epsilon for val in denom]
            hlist = list(hlist)

            denom = list(denom)
            loglist = np.log(2*epsilon + np.divide(hlist,denom))

            negl = np.dot(y,loglist)-np.sum(np.log(2*epsilon + denom))
            return negl

        def negjac(beta_est):
            hlist = h(beta_est,X)
            njac = -1*((np.subtract(y,hlist).transpose()@X))
            return njac


        beta_est = np.zeros(Ncov)
        #betatries = np.append(betatries)minimize(negloglikelihood,beta_est,method='Nelder-Mead')['x']
        #betatries[sample]=minimize(negloglikelihood,beta_est,method='Newton-CG', jac=negjac)['x']
        start_time = time.time()
        output = minimize(negloglikelihood,beta_est,method='Nelder-Mead')
        bTF = output['success']
        print(output)
        print('true beta is ', beta)
        print('ouput is', output['x'])
        print(bTF)
        elapsed_time = time.time()-start_time
        times.append([elapsed_time])
        if bTF == True:
            s=s+1

    probconv=s/N
    print(probconv, times)
    return probconv,betatries



#Gamma values to test
gamma = np.arange(0,10.4,0.4)

#Number of samples per square
N=1

#Number of observations per sample
n=400

#Distribution of gamma0/beta0. Set both to zero here. Change others' value appropriately in the loop.
beta0 = 0
gamma0 = 0

#Kappa value up to 0.6 inclusive and number of covariates:
kappa=np.arange(0,0.65,0.05)
p=n*kappa

#Array to store hMLE values
A=np.zeros((len(p),len(gamma)))
betadata={}
VN=scipy.stats.t.var(df)
countdown=len(kappa)*len(gamma)

#Obtain data
for j in range(len(gamma)):
    gamma0 = gamma[j]
    for i,Ncov in enumerate(p):
        print(countdown, gamma0, Ncov)
        countdown=countdown-1
        A[i,j], betatries = successprob(beta0,gamma0,Ncov,n,N,VN)
        kap=kappa[i]
#        bet = betatries
#        bet= pd.DataFrame(bet)
#        for i in range(np.shape(bet)[1]):
#            plt.plot(bet[:,i])
#            plt.show()



#Plot
pd.DataFrame(A).to_csv("array_NM_tn400_%d.csv" %N, header=None, index=None)
estimates=pd.DataFrame(betadata)
estimates.to_csv("betaestimates_NM_t400_%d.csv" %N)


plt.plot(betadata['kappa0.50_gamma4.0'])

#Save
plt.savefig('s')
...