Почему minimal_scalar не минимизирует правильно?

Question

Я новый пользователь Python, поэтому терпите меня, если этот вопрос очевиден.

Я пытаюсь найти значение lmbda, которое минимизирует следующую функцию, учитывая фиксированный вектор Z и скалярную сигма:

def sure_sft(z,lmbda, sigma):
     indicator = np.abs(z) <= lmbda;
     minimum = np.minimum(z**2,lmbda**2);
     return -sigma**2*np.sum(indicator) + np.sum(minimum);

Когда я передаю значения lmbda вручную, я обнаруживаю, что функция выдает правильное значение sure_stf.Однако, когда я пытаюсь использовать следующий код, чтобы найти значение lmbda, которое минимизирует sure_stf:

minimize_scalar(lambda lmbda: sure_sft(Z, lmbda, sigma))

, это дает мне неверное значение для sure_stf (-8.6731 для lmbda = 0.4916).Если я передам 0.4916 вручную в sure_sft, я получу -7.99809.Что я делаю неправильно?Буду признателен за любой совет!

РЕДАКТИРОВАТЬ: я вставил свой код ниже.Данные от: https://web.stanford.edu/~chadj/HallJones400.asc

import pandas as pd 
import numpy as np
from scipy.optimize import minimize_scalar

# FUNCTIONS

# Calculate orthogonal projection of g onto f  
def proj(f, g):  
    return  ( np.dot(f,g) / np.dot(f,f) ) * f   

def gs(X):

    # Copy of X -- will be used to store orthogonalization
    F = np.copy(X) 

    # Orthogonalize design matrix
    for i in range(1, X.shape[1]):            # Iterate over columns of X
        for j in range(i):                    # Iterate over columns less than current one
            F[:,i] -= proj(F[:,j], X[:,i])    # Subtract projection of x_i onto f_j for all j<i from F_i

    # normalize each column to have unit length
    norm_F=( (F**2).mean(axis=0) ) ** 0.5     # Row vector with sqrt root of average of the squares of each column
    W = F/norm_F                              # Normalize

    return W

# SURE for soft-thresholding    
def sure_sft(z,lmbda, sigma):
    indicator = np.abs(z) <= lmbda
    minimum = np.minimum(z**2,lmbda**2)
    return -sigma**2*np.sum(indicator) + np.sum(minimum)

# Import data.
data_raw =  pd.read_csv("hall_jones1999.csv")

# Drop missing observations.
data = data_raw.dropna(subset=['logYL', 'Latitude'])

Y = data['logYL']
Y = np.array(Y)
N = Y.size

# Create design matrix.
design = np.empty([data['Latitude'].size,15])

design[:,0] = 1
for j in range(1, 15):
    design[:,j] = data['Latitude']**j

K = design.shape[1]

# Use Gramm-Schmidt on design matrix.
W = gs(design)    
Z = np.dot(W.T, Y)/N    

# MLE
mu_mle = np.dot(W, Z)

# Soft-thresholding
# Use MLE residuals to calculate sigma for SURE calculation
sigma = np.sqrt(np.sum((Y - mu_mle)**2)/(N-K))

# Write SURE as a function of lmbda
sure = lambda lmbda: sure_sft(Z, lmbda, sigma)

# Find SURE-minimizing lmbda
lmbda = minimize_scalar(sure).x
min_sure = minimize_scalar(sure).fun #-8.673172212265738

# Compare to manually inputting minimized lambda into sure_sft
# I'm s
act_sure1 = sure_sft(Z, 0.49167598, sigma) #-7.998060514873529
act_sure2 = sure_sft(Z, 0.491675989, sigma) #-8.673172212306728

Answer 1

Вы на самом деле не делаете ничего плохого.Я только что проверил код и подтвердил, что lmbda имеет значение 0.4916759890416824 в конце скрипта.Вы можете убедиться в этом сами, добавив следующие строки в конец скрипта:

print(lmbda)
print(sure_sft(Z, lmbda, sigma))

при запуске скрипта вы должны увидеть:

0.4916759890416824
-8.673158394698172

Единственное, что яМожно предположить, что каким-то образом подпрограмма, которую вы использовали для распечатки lmbda, была настроена на печать только фиксированного числа цифр чисел с плавающей запятой, или каким-то образом распечатка была иным образом усечена.