Untitled

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the original function
def func(t, D, C, gamma):
    return D + C * (t ** gamma)

# Fixed parameters for synthetic data
D_true = 750
C_true = -573
gamma_true = 0.0324

# Let's define a function to calculate the root mean square error (RMSE) for the fit
def calculate_rmse(y_true, y_pred):
    return np.sqrt(np.mean((y_true - y_pred) ** 2))

best_rmse = float('inf')
best_D = None
best_params_log = None

# Generate synthetic data using more data points
t_values = np.linspace(50, 200, 1000)
y_values = func(t_values, D_true, C_true, gamma_true)
noise = np.random.normal(scale=0.01, size=len(t_values))
y_values_noise = y_values + noise

# Loop over possible D values to find the one that minimizes the RMSE
for D in np.linspace(749, 751, 100000):
    # Log-transform the data considering the 'D' term
    y_values_noise_log = np.ma.log(np.ma.masked_less_equal(-y_values_noise + D, 0))
    
    # Mask both t and x arrays to only consider valid log-transformed points
    common_mask = np.ma.mask_or(np.ma.getmask(t_values_log), np.ma.getmask(y_values_noise_log))
    t_values_log_masked = np.ma.masked_array(t_values_log, mask=common_mask).compressed()
    y_values_noise_log_masked = np.ma.masked_array(y_values_noise_log, mask=common_mask).compressed()

    # Fit the curve to the log-transformed data
    params_log, _ = curve_fit(func_log, t_values_log_masked, y_values_noise_log_masked, p0=[np.log(200), 0.5], maxfev=500000)

    # Generate curve for the log-transformed fitted data
    y_fit_log = func_log(t_values_log_masked, *params_log)

    # Calculate the RMSE for this fit
    rmse = calculate_rmse(y_values_noise_log_masked, y_fit_log)

    # Update the best fit if this one is better
    if rmse < best_rmse:
        best_rmse = rmse
        best_D = D
        best_params_log = params_log

# Extract the best fitted parameters
log_neg_C_fit, gamma_fit_log_transformed = best_params_log

# Recalculate C from log(-C)
C_fit_log_transformed = -np.exp(log_neg_C_fit)

best_D, C_fit_log_transformed, gamma_fit_log_transformed, best_rmse