[Python] Fitting with Maximum likelihood estimation in python returns initial parameters

I would like to use minimize from scipy.optimize to fit my data with a specific function:

I defined this log likelihood to this aim:

def log_likelihood(theta, freq, CPSD, CPSD_err):
'''
This function computes the log-likelihood for a given set of model parameters θ in a Bayesian framework,
comparing the observed data (freq, CPSD, CPSD_err) to a model.
- theta: A tuple containing the model parameters:
- freq: The independent variable eg frequency
- CPSD: The observed dependent variable: cross power spectral density
- CPSD_err: The uncertainties (errors) on the observed values of yy.
'''

#extract the model parameters
alpha_0, freq_O, alpha, CPSD_err = theta

# Define the model
model = alpha_0 * (freq_O**2 + freq**2)**(-(alpha + 3/2))

# Calculate the total variance:
# yerr2: The squared measurement uncertainties.
# exp(2⋅logf)=f2, which represents the fractional noise.
sigma2 = CPSD_err**2 #+ model**2 * np.exp(2 * log_f)

# computes and returns the log-likelihood:
# This is based on the formula for the log of a Gaussian likelihood:
return -0.5 * np.sum((CPSD - model) ** 2 / sigma2 + np.log(sigma2))

# step 2: Define the Negative Log-Likelihood
def nll(*args):
'''
Define the Negative Log-Likelihood
Minimizing the negative log-likelihood is equivalent to maximizing the log-likelihood.
Purpose: This function will be minimized to estimate the model parameters.
'''
return -log_likelihood(*args)


following this tutorial: https://emcee.readthedocs.io/en/stable/tutorials/line/

This is my function to do this:

def fit_CPSD_flattened(freq, CPSD, CPSD_err, freq_max, alpha_0, alpha, freq_O):
'''
This function fit the power spectrum with a flattened spectrum
- freq: frequencies, my data for the x axis
- CPSD: cross power spectral density, my data for the y axis
- freq-Max: freq until which I want to fit
- alpha_0: initial guess for parameter alpha_0
- alpha: initial guess for parameter alpha
- freq_0: initial guess for parameter freq_0
'''

# find the index of the maximum frequency for the fitting
# np.where returns an array with elements from x where condition is True, and elements from y elsewhere.
idx_high = np.where(freq >= freq_max)[0][0]
print(f'{freq[idx_high]} Hz')

# select the data range that will be fitted
CPSD_selected = CPSD[:idx_high]
freq_selected = freq[:idx_high]

CPSD_err_scalar = np.std(CPSD_err)

# initial guess for the parameters
initial = np.array([alpha_0, freq_O, alpha, CPSD_err_scalar])


# Define bounds for the parameters
# Each tuple corresponds to (min, max) for a parameter
bounds = [
(1e-11, 1e-4), # Bounds for alpha_0
(0.0001, 2.0), # Bounds for freq_O
(-6, 0), # Bounds for alpha
(-np.inf, np.inf) # Bounds for CPSD_err
]

# step 4: Minimize the Negative Log-Likelihood
soln = minimize(nll, initial, args=(freq_selected, CPSD_selected, CPSD_err_scalar), method='L-BFGS-B', bounds=bounds)

# Step 5: Extract the Best-Fit Parameters
alpha_0_ml, freq_O_ml, alpha_ml, CPSD_err_ml = soln.x

return alpha_0_ml, freq_O_ml, alpha_ml


and this is my code where I call the function with my data:

# std of the CPSD std
CPSD_err_scalar = np.std(CPSD_std_19)

#initial guess for the parameter
alpha_0 = 5e-10
freq_O = 0.006
alpha = -0.30

freq_max = 1 #Hz

# In this specific example:
alpha_0_ml, freq_O_ml, alpha_ml = fit_CPSD_flattened(bins_19[CPSD_mean_19>0], CPSD_mean_19[CPSD_mean_19>0], CPSD_err_scalar, freq_max, alpha_0, alpha, freq_O)

print(f'alpha_0_ml = {alpha_0_ml}')
print(f'freq_O_ml = {freq_O_ml}')
print(f'alpha_ml = {alpha_ml}')


slope = alpha_ml + 3/2

plt.plot(bins_19[CPSD_mean_19>0], CPSD_mean_19[CPSD_mean_19>0])
# fit
plt.plot(bins_19, flattening_CPSD(bins_19, alpha_0_ml, freq_O_ml, alpha_ml), color='lime',
label=f'Fit: 1/k_out = {np.round(freq_O_ml, 4)}Hz, slope = -({np.round(slope*2, 3)})')
# label=f'Fit: 1/k_out = ({np.round(freq_O_ml, 4)} ± {np.round(freq_O_ml_err, 4)})Hz, slope = -({np.round(slope*2, 3)} ± {np.round(slope_err, 3)})')
plt.loglog()
plt.legend()


but this ONLY returns the initial parameter... In the picture, the green fit is completely off the data in blue, and the fit only takes the value I gave for the initial parameter; this is obviously wrong. This drives me crazy, I don't know how to use this scipy minimize correctly, could someone help me?



This is the whole code for a reproducible example:

import numpy as np
from numpy import loadtxt
import matplotlib.pyplot as plt
from astropy.io import fits
import datetime
import scipy.stats as stats
from scipy import signal
from scipy.fftpack import fft, fftfreq
import math
from astropy.time import Time #convert julian date to normal time
from scipy.signal import find_peaks # try the find peaks scipy function
import healpy as hp # to plot the QUIJOTE's map
import matplotlib.ticker as ticke
from scipy.optimize import curve_fit # to fit functions
from matplotlib import cm # for color map
import os # to create directories
from scipy.interpolate import UnivariateSpline
import pandas as pd
import sys
from datetime import datetime
from astropy.coordinates import EarthLocation, AltAz, SkyCoord
from astropy import units as u
from tqdm import tqdm # print a progression bar
import time # use a timer
from datetime import datetime, timezone
from matplotlib.ticker import MultipleLocator # to sert the position of minor ticks
import matplotlib as mpl # to increase the thickness of plot frame
from scipy.optimize import minimize


def PS_oof_fct02(freq, f_s, signoise, alpha, f_knee):
'''
Generate the theoretical power spectrum of a 1/f noise signal.
Same as above just I deleted the one in PS of the 1/f signal.
- freq: frequency range of the power spectrum
- f_s: sampling frequency of the theoretical signal
- signoise: level of the white noise signal: fixed!
- alpha: slope of the 1/f noise signal
_ f_knee: knee frequency of the 1/f noise signal
'''
############################# 1) Power spectrum of the white noise
# number of sample in my theoretical observation
nsample = len(freq)

# generate white noise level white std signoise
White_noise = np.random.normal(0, signoise, nsample) # loc, scale, size

# Computing the one-dimensional discrete Fourier Transform of the white noise. Normalised
White_noise_fft = direct_fft(White_noise, f_s) #direct_fft(f,f_s): np.fft.fft(White_noise) / f_sampling

# Power spectrum of the white noise signal
PS_white_noise = np.real(White_noise_fft * np.conjugate(White_noise_fft))

############################# 2) Creates the correction term for the white noise which will be
# the 1/f noise
# array of zeros with the same size as F to be filled with the 1/f correction term
PS_corr = freq * 0

############################# 3) fill the array for the positiv frequencies
#mask to select only the positiv frequency
mask_pos = (freq > 0)

# apply the mask to the array of frequencies
freq_pos = freq[mask_pos]

# define the corrrection term for the 1/f noise for positive frequencies
PS_corr[mask_pos] = (f_knee / freq_pos) ** alpha

############################# 4) fill the array for the negativ frequencies
#mask to select only the negative frequencies
mask_neg = (freq < 0)

# apply the mask to the array of frequencies
freq_neg = freq[mask_neg]

#define the corrrection term for the 1/f noise for negative frequencies
PS_corr[mask_neg] = (f_knee / np.abs(freq_neg)) ** alpha

############################# 5) fill the array for the frequencies equal zero
# mask to select frequencies equal to zeros
mask_zero = (freq == 0)

#define the corrrection term for the 1/f noise for frequencies equal to zero
PS_corr[mask_zero] = (f_knee / np.min(freq[freq > 0])) ** alpha

############################# 6) Power spectrum of the one over f noise signal
# multiply the white noise level with the correction term
PS_oof = PS_white_noise * PS_corr

############################# 7) 1/f noise signal in time domain
# fft of the 1/f noise signal to get the time domain signal
oof_fft = White_noise_fft * np.sqrt(PS_corr)

# then do in inverse FFT of it to obtain it in time domain
oof_time = np.fft.ifft(oof_fft).real * f_s

return PS_oof, oof_time, PS_corr, White_noise_fft, White_noise

def PS_white_noise_fct(freq, f_s, signoise):
'''
Generate the theoretical power spectrum of a white noise signal
- freq: frequency range of the power spectrum
- f_s: sampling frequency of the theoretical signal
- signoise: amplitude of the white noise level
'''

# number of sample in my theoretical observation
nsample = len(freq)

# generate white noise level white std signoise
White_noise_time = np.random.normal(0, signoise, nsample) # loc, scale, size

# Computing the one-dimensional discrete Fourier Transform of the white noise. Normalised
White_noise_fft = direct_fft(White_noise_time, f_s) #direct_fft(f,f_s): np.fft.fft(White_noise) / f_sampling

# Power spectrum of the white noise signal
PS_white_noise = np.real(White_noise_fft * np.conjugate(White_noise_fft))

return PS_white_noise, White_noise_time

def direct_fft(f,f_s):
return np.fft.fft(f) / f_s

def log_likelihood(theta, freq, CPSD, CPSD_err):
'''
This function computes the log-likelihood for a given set of model parameters θ in a Bayesian framework,
comparing the observed data (freq, CPSD, CPSD_err) to a model.
- theta: A tuple containing the model parameters:
- freq: The independent variable eg frequency
- CPSD: The observed dependent variable: cross power spectral density
- CPSD_err: The uncertainties (errors) on the observed values of yy.
'''

#extract the model parameters
alpha_0, freq_O, alpha, CPSD_err = theta

# Define the model
model = alpha_0 * (freq_O**2 + freq**2)**(-(alpha + 3/2))

# Calculate the total variance:
# yerr2: The squared measurement uncertainties.
# exp(2⋅logf)=f2, which represents the fractional noise.
sigma2 = CPSD_err**2 #+ model**2 * np.exp(2 * log_f)

# computes and returns the log-likelihood:
# This is based on the formula for the log of a Gaussian likelihood:
return -0.5 * np.sum((CPSD - model) ** 2 / sigma2 + np.log(sigma2))

# step 2: Define the Negative Log-Likelihood
def nll(*args):
'''
Define the Negative Log-Likelihood
Minimizing the negative log-likelihood is equivalent to maximizing the log-likelihood.
Purpose: This function will be minimized to estimate the model parameters.
'''
return -log_likelihood(*args)


def fit_CPSD_flattened(freq, CPSD, CPSD_err, freq_max, alpha_0, alpha, freq_O):
'''
This function fit the power spectrum with a flattened spectrum
- freq: frequencies, my data for the x axis
- CPSD: cross power spectral density, my data for the y axis
- freq-Max: freq until which I want to fit
- alpha_0: initial guess for parameter alpha_0
- alpha: initial guess for parameter alpha
- freq_0: initial guess for parameter freq_0
'''

# find the index of the maximum frequency for the fitting
# np.where returns an array with elements from x where condition is True, and elements from y elsewhere.
idx_high = np.where(freq >= freq_max)[0][0]
print(f'{freq[idx_high]} Hz')

# select the data range that will be fitted
CPSD_selected = CPSD[:idx_high]
freq_selected = freq[:idx_high]

CPSD_err_scalar = np.std(CPSD_err)

# initial guess for the parameters
initial = np.array([alpha_0, freq_O, alpha, CPSD_err_scalar])


# Define bounds for the parameters
# Each tuple corresponds to (min, max) for a parameter
bounds = [
(1e-11, 1e-4), # Bounds for alpha_0
(0.0001, 2.0), # Bounds for freq_O
(-6, 0), # Bounds for alpha
(-np.inf, np.inf) # Bounds for CPSD_err
]

# step 4: Minimize the Negative Log-Likelihood
soln = minimize(nll, initial, args=(freq_selected, CPSD_selected, CPSD_err_scalar), method='L-BFGS-B', bounds=bounds)

# Step 5: Extract the Best-Fit Parameters
alpha_0_ml, freq_O_ml, alpha_ml, CPSD_err_ml = soln.x

return alpha_0_ml, freq_O_ml, alpha_ml


def flattening_CPSD(freq, alpha_0, freq_O, alpha):
'''
Function I want to fit to my data
- freq: x axis
- alpha_0, freq_O, alpha: functions parameters
'''
return alpha_0 * (freq_O**2 + freq**2)**(-(alpha + 3/2))

############################################################## This part generate synthetic data
# sampling frequency
f_sampling_MFI = 1/0.04 # 25 Hz

nsamp = 4000

# standard deviation for the repetition of the 1/f noise
signoise = 0.01

# standard deviation for the repetition of the white noise signal 1
signoise01 = 0.3

# standard deviation for the repetition of the white noise signal 2
signoise02 = 0.3
alpha = 2.66
f_knee = 5.0
f_s = f_sampling_MFI

# define the frequencies for this sample size ans spacing
freq = np.fft.fftfreq(nsamp, d = 1/f_sampling_MFI)

# theoretical atmospheric 1/f noise
oof_ps, time_oof, ps_corr, white_noise_fft, white_noise = PS_oof_fct02(freq, f_sampling_MFI, signoise, alpha, f_knee)

# theoretical intrumental noise horn 2
PS_white_noise01, time_white_noise01 = PS_white_noise_fct(freq, f_sampling_MFI, signoise01)

# Power spectrum of horn 2 theoretical white noise + atmopsheric 1/f noise level
PS_h2_theo = oof_ps + PS_white_noise01



# Shift the zero-frequency component to the center of the spectrum.
freq_shf = np.fft.fftshift(freq)
PS_h2_theo_shf = np.fft.fftshift(PS_h2_theo)


plt.figure()
# PS horn 2 theo
plt.plot(freq_shf, PS_h2_theo_shf)
plt.loglog()
plt.xlabel('Freq [Hz]')
plt.ylabel('Power')


############################################################## This part tries to fit the data

# std of the CPSD std
CPSD_err_scalar = np.std(PS_h2_theo)

#initial guess for the parameter
alpha_0 = 5e-10
freq_O = 0.006
alpha = -0.30

freq_max = 1 #Hz

# In this specific example:
alpha_0_ml, freq_O_ml, alpha_ml = fit_CPSD_flattened(freq_shf, PS_h2_theo_shf, CPSD_err_scalar, freq_max, alpha_0, alpha, freq_O)

print(f'alpha_0_ml = {alpha_0_ml}')
print(f'freq_O_ml = {freq_O_ml}')
print(f'alpha_ml = {alpha_ml}')

# parameter to print in the figure
slope = alpha_ml + 3/2

plt.figure()
plt.plot(freq_shf, PS_h2_theo_shf)
# fit
plt.plot(freq_shf, flattening_CPSD(freq_shf, alpha_0_ml, freq_O_ml, alpha_ml), color='lime',
label=f'Fit: 1/k_out = {np.round(freq_O_ml, 4)}Hz, slope = -({np.round(slope*2, 3)})')
plt.loglog()
plt.legend()

Continue reading...