OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

How to choose the best fitting method?

  • Thread starter Thread starter sam wolfe
  • Start date Start date
S

sam wolfe

Guest
I want to fit parameters x to a function fp, given by

enter image description here

where v and L are given parameters. Computing fp can be done relatively quickly by vectorizing it and using some uniform filtering. The following code implements this optimization problem, where the function cfit (through fitf) slowly adjusts the values of xs by calculating fp at each step and comparing it to the real data list, according to

enter image description here

Question: Is there a way of improving this fitting method for this specific function, and the way I am computing it? I tried playing around with scipy.optimize, but everything seems slower than this simple multiplicative gradient method.

Here is the code implementation, where the mean-square error is tracked

Code:
import numpy as np
import math

def fitfunction(list, maxiter, error_threshold):
    
    v = 1.4
    st = 5
    exp_v = np.exp(-1/v)
    x00 = np.array([(math.pi/(4*v))*i**(-2) for i in list])
    lm = 1000
    
    def mse(y_true, y_pred):
        mse_value = sum((yt - yp) ** 2 for yt, yp in zip(y_true, y_pred)) / len(y_true)
        return mse_value
    
    def fast_roll_add(dst, src, shift):
        dst[shift:] += src[:-shift]
        dst[:shift] += src[-shift:]
    
    # Main function
    def fp(x, L, v):
        n = len(x)
        y = np.zeros(n)
        last_exp_2_raw = np.zeros(n)
        last_exp_2 = np.ones(n)
        unitary = x.copy()
        for k in range(L+1):
            if k != 0:
                fast_roll_add(unitary, x, k)
                fast_roll_add(unitary, x, -k)
            exp_1_raw = last_exp_2_raw
            exp_1 = last_exp_2
            exp_2_raw = exp_1_raw + unitary / v
            exp_2 = np.exp(-exp_2_raw)
            y += (exp_1 - exp_2) / unitary
            last_exp_2_raw = exp_2_raw
            last_exp_2 = exp_2
        return y

    # Fitting functions
    def fitf(time, lst, x0, j):
        return x0[j] * (lst[j] / time[j])**2
    def cfit(time, lst, x0):
        result = np.empty_like(x0)
        for j in range(len(x0)):
            if fitf(time, lst, x0, j) < 10**(-20):
                result[j] = 10**(-20)
            elif abs(time[j] - lst[j]) < .5:
                result[j] = x0[j]
            else:
                result[j] = fitf(time, lst, x0, j)
        return result
    
    xs = x00
    ys = fp(xs, len(xs)//st, v)
    err = 10**10
    
    for j in range(maxiter):
        if err > error_threshold:
            xs = cfit(list, ys, xs) # Fitting
            ys = fp(xs, len(xs)//st, v) # Update function values
            err = mse(list[lm:-lm], ys[lm:-lm])
            print(str(j+1) + '/' + str(maxiter) + ' err: ' + str('{:.20f}'.format(err)), end="\r")            
        else:
            break

    fire_rates = ['{:.20f}'.format(i) for i in xs]
    time_sim = ys
    
    return [fire_rates, time_sim]

And some working example, with generated data

Code:
from scipy.ndimage import gaussian_filter1d

# Data generation
np.random.seed(18)
data = np.random.normal(loc=0, scale=10, size=10000).cumsum()
data = (data - data.min()) / (data.max() - data.min()) * 500

# Gaussian smoothing
data = gaussian_filter1d(data, sigma=100)

data_fit = fitfunction(data, 100, 2)

and to visualize it

Code:
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.plot(data, label='Data')
plt.plot(data_fit[1], label='Data Fit')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Data and Data Fit')
plt.legend()
plt.show()

enter image description here

Brief comments:

  • Naturally, this fit "seems" quite good and relatively fast, but I plan to do it on a larger scale (10^6 points), and I also want to guarantee the best possible fit. My multiplicative method might be limiting.
  • The implementation of fp is periodic in nature, which leads to massive errors in the ends. Hence, the error err is calculated excluding these regions.
  • The fitting function provides an even better fit

Code:
def fitf(time, lst, x0, j):
    return x0[j]**(np.log(time[j]) / np.log(lst[j]))

However, when the values come too close to the real value, it seems to be highly unstable.
<p>I want to fit parameters <code>x</code> to a function <code>fp</code>, given by</p>
<p><a href="https://i.sstatic.net/jyCwdRBF.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jyCwdRBF.png" alt="enter image description here" /></a></p>
<p>where <code>v</code> and <code>L</code> are given parameters. Computing <code>fp</code> can be done relatively quickly by vectorizing it and using some uniform filtering. The following code implements this optimization problem, where the function <code>cfit</code> (through <code>fitf</code>) slowly adjusts the values of <code>xs</code> by calculating <code>fp</code> at each step and comparing it to the real data <code>list</code>, according to</p>
<p><a href="https://i.sstatic.net/Tnya1WJj.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Tnya1WJj.png" alt="enter image description here" /></a></p>
<p><strong>Question:</strong> Is there a way of improving this fitting method for this specific function, and the way I am computing it? I tried playing around with <code>scipy.optimize</code>, but everything seems slower than this simple multiplicative gradient method.</p>
<p>Here is the code implementation, where the mean-square error is tracked</p>
<pre><code>import numpy as np
import math

def fitfunction(list, maxiter, error_threshold):

v = 1.4
st = 5
exp_v = np.exp(-1/v)
x00 = np.array([(math.pi/(4*v))*i**(-2) for i in list])
lm = 1000

def mse(y_true, y_pred):
mse_value = sum((yt - yp) ** 2 for yt, yp in zip(y_true, y_pred)) / len(y_true)
return mse_value

def fast_roll_add(dst, src, shift):
dst[shift:] += src[:-shift]
dst[:shift] += src[-shift:]

# Main function
def fp(x, L, v):
n = len(x)
y = np.zeros(n)
last_exp_2_raw = np.zeros(n)
last_exp_2 = np.ones(n)
unitary = x.copy()
for k in range(L+1):
if k != 0:
fast_roll_add(unitary, x, k)
fast_roll_add(unitary, x, -k)
exp_1_raw = last_exp_2_raw
exp_1 = last_exp_2
exp_2_raw = exp_1_raw + unitary / v
exp_2 = np.exp(-exp_2_raw)
y += (exp_1 - exp_2) / unitary
last_exp_2_raw = exp_2_raw
last_exp_2 = exp_2
return y

# Fitting functions
def fitf(time, lst, x0, j):
return x0[j] * (lst[j] / time[j])**2
def cfit(time, lst, x0):
result = np.empty_like(x0)
for j in range(len(x0)):
if fitf(time, lst, x0, j) < 10**(-20):
result[j] = 10**(-20)
elif abs(time[j] - lst[j]) < .5:
result[j] = x0[j]
else:
result[j] = fitf(time, lst, x0, j)
return result

xs = x00
ys = fp(xs, len(xs)//st, v)
err = 10**10

for j in range(maxiter):
if err > error_threshold:
xs = cfit(list, ys, xs) # Fitting
ys = fp(xs, len(xs)//st, v) # Update function values
err = mse(list[lm:-lm], ys[lm:-lm])
print(str(j+1) + '/' + str(maxiter) + ' err: ' + str('{:.20f}'.format(err)), end="\r")
else:
break

fire_rates = ['{:.20f}'.format(i) for i in xs]
time_sim = ys

return [fire_rates, time_sim]
</code></pre>
<p>And some working example, with generated data</p>
<pre><code>from scipy.ndimage import gaussian_filter1d

# Data generation
np.random.seed(18)
data = np.random.normal(loc=0, scale=10, size=10000).cumsum()
data = (data - data.min()) / (data.max() - data.min()) * 500

# Gaussian smoothing
data = gaussian_filter1d(data, sigma=100)

data_fit = fitfunction(data, 100, 2)
</code></pre>
<p>and to visualize it</p>
<pre><code>import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.plot(data, label='Data')
plt.plot(data_fit[1], label='Data Fit')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Data and Data Fit')
plt.legend()
plt.show()
</code></pre>
<p><a href="https://i.sstatic.net/zOboLGx5.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zOboLGx5.png" alt="enter image description here" /></a></p>
<p><strong>Brief comments:</strong></p>
<ul>
<li>Naturally, this fit "seems" quite good and relatively fast, but I plan to do it on a larger scale (10^6 points), and I also want to guarantee the best possible fit. My multiplicative method might be limiting.</li>
<li>The implementation of <code>fp</code> is periodic in nature, which leads to massive errors in the ends. Hence, the error <code>err</code> is calculated excluding these regions.</li>
<li>The fitting function provides an even better fit</li>
</ul>
<pre><code>def fitf(time, lst, x0, j):
return x0[j]**(np.log(time[j]) / np.log(lst[j]))
</code></pre>
<p>However, when the values come too close to the real value, it seems to be highly unstable.</p>
 

Latest posts

Top