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

问题描述

您好，如标题所示，我一直在尝试为PDF添加适合的指数和幂定律。 如下图所示：

我正在使用的代码生成底层图形：

代码是这样的：

   a11=[9.76032106e-02,6.73754187e-02,3.20683249e-02,2.21788509e-02,2.70850237e-02,9.90377323e-03,2.11573411e-02,8.46232347e-03,8.49027869e-03,7.33997745e-03,5.71819070e-03,4.62720448e-03,4.11562884e-03,3.20064313e-03,2.66192941e-03,1.69116510e-03,1.94355212e-03,2.55224949e-03,1.23822395e-03,5.29618250e-04,4.03769641e-04,3.96865740e-04,3.38530868e-04,2.04124701e-04,1.63913557e-04,2.04486864e-04,1.82216592e-04,1.34708400e-04,9.24289261e-05,9.55074181e-05,8.13695322e-05,5.15610541e-05,4.15425149e-05,4.68101099e-05,3.33696885e-05,1.61893058e-05,9.61743970e-06,1.17314090e-05,6.65239507e-06]

b11=[3.97213201e+00,4.77600082e+00,5.74255432e+00,6.90471618e+00,8.30207306e+00,9.98222306e+00,1.20023970e+01,1.44314081e+01,1.73519956e+01,2.08636432e+01,2.50859682e+01,3.01627952e+01,3.62670562e+01,4.36066802e+01,5.24316764e+01,6.30426504e+01,7.58010432e+01,9.11414433e+01,1.09586390e+02,1.31764173e+02,1.58430233e+02,1.90492894e+02,2.29044305e+02,2.75397642e+02,3.31131836e+02,3.98145358e+02,4.78720886e+02,5.75603061e+02,6.92091976e+02,8.32155588e+02,1.00056488e+03,1.20305636e+03,1.44652749e+03,1.73927162e+03,2.09126048e+03,2.51448384e+03,3.02335795e+03,3.63521656e+03,4.37090138e+03]
                                                      
    plt.plot(b11,a11,'ro')
    plt.yscale("log")
    plt.xscale("log")
    
    plt.show()
     

我想基于卡方误差最小化方法，在更小的时间上将幂定律拟合到下图上，并在更长的时间上添加指数拟合。

以csv格式保存的x轴数据：

x轴的数据：

解决方法

正如我在评论中所提到的，我认为您可以通过一个常数项将幂律和指数相结合。另外，数据看起来可以由两个幂定律拟合。尽管评论表明确实存在指数行为。无论如何，我在这里展示两种方法。在这两种情况下，我都尽量避免使用任何类型的分段定义。这样还可以确保$ C ^ infty $。

在第一种方法中，a * x**( -b )代表小x，而a1 * exp( -d * x )代表大x。想法是选择一个c，使得对于所需的小c，幂律比x大得多，但对于其他情况，则要小得多。 这允许使用我的注释中提到的功能，即( a * x**( -b ) + c ) * exp( -d * x ) 。可以将c作为过渡参数。

在替代方法中，我采用了两个幂律。因此，存在两个区域，在第一个功能中，一个区域较小，在第二个区域中，第二个区域较小。因为我一直想要较小的函数，所以我进行了求反，即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,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,1e7,4,1 ]
atesty = np.fromiter( ( 10**double_power_log( x,*aguess ) for x in testx ),np.float )
asol,_ = curve_fit( double_power_log,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,label="guess p-p" )
ax.plot( testx,fity,ls='-',label="fit p-e: {}".format( sol ) )
ax.plot( testx,afity,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() 

结果看起来像

,

为了完整起见，我想添加一个具有分段定义的解决方案。因为我希望函数连续且可微，所以指数定律的参数不是完全自由的。使用f = a * x**(-b)和g = alpha * exp( -beta * x )并在x0进行过渡时，我选择( a,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",' )

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

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

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

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


def alpha_beta( x,x0 ):
    """
    continuous and differentiable define alpha and beta
    free parameter is the point where I connect
    """
    sol = minimize(alpha_fun,.005,args=( a,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,b ) )/ x0 )
    ### 
    if isinstance( x,( np.ndarray,list,tuple ) ):
        out = list()
        for v in x:
            out.append( p_w( v,x0 ) )
    else:
        out = p_w( v,x0 )
    return out

sol,_ = curve_fit( alpha_beta,p0=[ .1,70. ] )

alpha0 = minimize(alpha_fun,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,100)
yl = alpha_beta( xl,*sol )

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

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

最终提供