Vfrdtyky

I am a new Python user, so bear with me if this question is obvious.

I am trying to find the value of lmbda that minimizes the following function, given a fixed vector Z and scalar sigma:

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);

When I pass in values of lmbda manually, I find that the function produces the correct value of sure_stf. However, when I try to use the following code to find the value of lmbda that minimizes sure_stf:

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

it gives me an incorrect value for sure_stf (-8.6731 for lmbda = 0.4916). If I pass in 0.4916 manually to sure_sft, I obtain -7.99809 instead. What am I doing incorrectly? I would appreciate any advice!

EDIT: I've pasted my code below. The data is from: 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

You're actually not doing anything wrong. I just tested out the code and confirmed that lmbda has a value of 0.4916759890416824 at the end of the script. You can confirm this for yourself by adding the following lines to the bottom of your script:

print(lmbda)

print(sure_sft(Z, lmbda, sigma))

when you run your script you should then see:

0.4916759890416824

-8.673158394698172

The only thing I can figure is that somehow the routine you were using to print out lmbda was set up to only print a fixed number of digits of floating point numbers, or somehow the printout was otherwise truncated.