在我的PDF中添加基于卡方误差最小化的幂律和指数拟合

2条回答

网友

1楼 · 编辑于 2024-05-14 10:01:20

正如我在评论中提到的，我认为你可以通过一个常数项来耦合幂律和指数。或者，数据看起来可以用两个幂律拟合。尽管这些评论表明确实存在指数行为。无论如何，我在这里展示了这两种方法。在这两种情况下，我都试图避免任何类型的分段定义。这也确保了$C^infty$

在第一种方法中，我们用a * x**( -b )表示小x，用a1 * exp( -d * x )表示大x。我们的想法是选择一个c，这样，对于所需的小x，幂律比c大得多，但在其他方面要小得多。这允许我的评论中提到的函数，即( a * x**( -b ) + c ) * exp( -d * x ) 。可以考虑{{CD5}}作为转换参数。

在替代方法中，我采用两个幂律。因此，有两个区域，第一个区域较小，第二个区域较小。因为我总是想要较小的函数，所以我进行逆求和，即f = 1 / ( 1 / f1 + 1 / f2 )。从下面的代码中可以看出，我添加了一个额外的参数（技术上是在[0，infty[）。该参数控制过渡的平滑度

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

data = np.loadtxt( "7jyRi.txt", delimiter=',' )

#### p-e: power and exponential coupled via a small constant term
def func_log( x, a, b, c, d ):
    return np.log10( ( a * x**( -b ) + c ) * np.exp( -d * x ) )

guess = [.1, .8, 0.01, .005 ]
testx = np.logspace( 0, 3, 150 )
testy = np.fromiter( ( 10**func_log( x, *guess ) for x in testx ), np.float )
sol, _ = curve_fit( func_log, data[ ::, 0 ], np.log10( data[::,1] ), p0=guess )
fity = np.fromiter( ( 10**func_log( x, *sol ) for x in testx ), np.float )

#### p-p: alternatively using two power laws
def double_power_log( x, a, b, c, d, k ):
    s1 = ( a * x**( -b ) )**k
    s2 = ( c * x**( -d ) )**k
    out = 1.0 / ( 1.0 / s1 + 1.0 / s2 )**( 1.0 / k )
    return np.log10( out )

aguess = [.1, .8, 1e7, 4, 1 ]
atesty = np.fromiter( ( 10**double_power_log( x, *aguess ) for x in testx ), np.float )
asol, _ = curve_fit( double_power_log, data[ ::, 0 ], np.log10( data[ ::, 1 ] ), p0=aguess )
afity = np.fromiter( ( 10**double_power_log( x, *asol ) for x in testx ), np.float )

#### plotting
fig = mp.figure( figsize=( 10, 8 ) )
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( data[::,0], data[::,1] ,ls='', marker='o', label="data" )
ax.plot( testx, testy ,ls=':', label="guess p-e" )
ax.plot( testx, atesty ,ls=':',label="guess p-p" )
ax.plot( testx, fity ,ls='-',label="fit p-e: {}".format( sol ) )
ax.plot( testx, afity ,ls='-', label="fit p-p: {}".format( asol ) )
ax.set_xscale( "log" )
ax.set_yscale( "log" )
ax.set_xlim( [ 5e-1, 2e3 ] )
ax.set_ylim( [ 1e-5, 2e-1 ] )
ax.legend( loc=0 )
mp.show()

结果看起来像

网友

2楼 · 编辑于 2024-05-14 10:01:20

为了完整性，我想添加一个具有分段定义的解决方案。因为我想要函数是连续的和可微的，指数定律的参数不是完全自由的。使用f = a * x**(-b)和g = alpha * exp( -beta * x )以及x0处的转换，我选择( a, b, x0 )作为自由参数。从这个alpha和beta开始。然而，这些方程没有简单的解，因此这本身就需要最小化

import matplotlib.pyplot as mp
import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import minimize
from scipy.special import lambertw

data = np.loadtxt( "7jyRi.txt", delimiter=',' )

def pwl( x, a, b):
    return a * x**( -b )

def expl( x, a, b ):
    return a * np.exp( -b * x )

def alpha_fun(alpha, a, b, x0):
    out = alpha - pwl( x0, a, b ) * expl(1, 1, lambertw( pwl( x0, -a * b/ alpha, b ) ) )
    return 1e10 * np.abs( out )**2

def p_w( v, a,b, alpha, beta, x0 ):
    if v < x0:
        out = pwl( v, a, b )
    else:
        out = expl( v, alpha, beta )
    return np.log10( out )


def alpha_beta( x, a, b, x0 ):
    """
    continuous and differentiable define alpha and beta
    free parameter is the point where I connect
    """
    sol = minimize(alpha_fun, .005, args=( a, b, x0 ) )### attention, strongly depends on starting guess, i.e might be a catastrophic fail
    alpha = sol.x[0]
    # ~print alpha
    beta = np.real( -lambertw( pwl( x0, -a * b/ alpha, b ) )/ x0 )
    ### 
    if isinstance( x, ( np.ndarray, list, tuple ) ):
        out = list()
        for v in x:
            out.append( p_w( v, a, b, alpha, beta, x0 ) )
    else:
        out = p_w( v, a, b, alpha, beta, x0 )
    return out

sol,_ = curve_fit( alpha_beta, data[ ::, 0 ], np.log10( data[ ::, 1 ] ), p0=[ .1, .8, 70. ] )

alpha0 = minimize(alpha_fun, .005, args=tuple(sol ) ).x[0]
beta0 = np.real( -lambertw( pwl( sol[2], -sol[0] * sol[1]/ alpha0, sol[1] ) )/ sol[2] )

xl = np.logspace(0,3,100)
yl = alpha_beta( xl, *sol )

pl = pwl( xl, sol[0], sol[1] )
el = expl( xl, alpha0, beta0 )
#### plotting
fig = mp.figure( figsize=( 10, 8 ) )
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( data[::,0], data[::,1] ,ls='', marker='o', label="data" )
ax.plot( xl, pl ,ls=':', label="p" )
ax.plot( xl, el ,ls=':', label="{:0.3e} exp(-{:0.3e} x)".format(alpha0, beta0) )
ax.plot( xl, [10**y for y in yl] ,ls='-', label="sol: {}".format(sol) )

ax.axvline(sol[-1], color='k', ls=':')
ax.set_xscale( "log" )
ax.set_yscale( "log" )
ax.set_xlim( [ 5e-1, 2e3 ] )
ax.set_ylim( [ 1e-5, 2e-1 ] )
ax.legend( loc=0 )
mp.show()

最终提供

相关问题更多 >

编程相关推荐

热门问题

热门文章