背景: 由于光时轨道效应，脉动变星的 O-C 图会有类似正弦的变化，在拟合过程中，会用考虑到轨道倾角的因素(sin i) 问题: 如何加入sin i 这一项来拟合呢？ 另: Oscillation modes: by a standard Fourier analysis of light curve . method: Period04 The main oscillation mode is ... and the second mode is ...., then the two oscillation mode were used to fitted the light curve . Reference: Jafarzadeh, S. J. 2017 New Astronomy. 2019. 07. 30
样条函数插值拟合 2014–02–11 09:26:49 在拟合势能函数的时候, 除解析式外, 也可以利用样条函数进行拟合. 样条拟合与其插值正好相反: 已知函数在节点上的值求任意位置的值, 做插值; 已知函数的某些组合值求函数节点上的值, 属于拟合. 由于样条函数可以化为节点函数值的线性表达式, 这样就可以将待求参数线性化, 得到最优情况下函数的形状, 为寻找合适的解析式提供依据, 当然也可以直接利用得到的离散数据拟合解析式. 样条函数可以是零阶, 一阶, 二阶, 三阶或更高阶. 实际使用中, 三阶使用最为普遍. 由于三次样条的构造需要求解一个三对角线性方程组, 其显式解很难得到, 所以线性化结果很繁琐. 零阶 在每一区间上样条函数为常量, 函数整体呈台阶状. 对等距情况, 计算时最好使用就近原则, 取最近点的值作为拟合点, 可用 i=nint(x/dx) 实现. 一阶 在每一区间上样条函数为线性函数, 函数整体呈折线状$.$ f i ( x ) = y i + y i + 1 − y i Δ x ( x − x i ) 此式自动满足函数值连续条件, 即零阶连续. 二阶 在每一区间上, 样条函数为二次函数, 整体一阶连续, 即有连续的导数. 但仅有节点的函数值不能唯一确定整个函数, 还须提供某一节点上的导数值, 一般可令端点的导数值为零. 二次样条函数在偶数点的曲率不连续. 由二阶开始, 插值函数不再具有局域性, 改变某一节点, 函数整体都会改变. 使得线性系数分离很困难. f i ( x ) k i + 1 = y i + k i ( x − x i ) + k i + 1 − k i 2 Δ x ( x − x i ) 2 = − k i + 2 y i + 1 − y i Δ x , k i = 2 y i + 1 − y i Δ x − k i + 1 可化简得 f i ( x ) α = y i + k i ( x − x i ) + ( y i + 1 − y i − k i Δ x ) ( x − x i Δ x ) 2 = α 2 y i + 1 + ( 1 − α 2 ) y i + ( 1 − α ) ω k i = ω / Δ x , ω = x − x i 对势能函数, 一般满足远距离处导数为零, 故可使用自然条件 k n = 0 , 由此, 可推知所有系数 k i . 令 k i = 2 Δ x T i , Δ i = y i + 1 − y i , 则 T i 满足递推式 T i = Δ i − T i + 1 可求得 T i = ∑ j = i n − 1 ( − 1 ) j − i Δ j 样条函数可写为 f i ( x ) = α 2 y i + 1 + ( 1 − α 2 ) y i + 2 α ( 1 − α ) T i , α = ( x − x i ) / Δ x 对不等距划分, 令 Δ i = y i + 1 − y i x i + 1 − x i , k i 满足如下递推式 k i = 2 Δ i − k i + 1 求得 k i = 2 ∑ j = i n − 1 ( − 1 ) j − i Δ j 三阶 对等距划分的均匀样条, 设节点为 1 , 2 , . . . . n , 若 x ∈ , a = x − x i , b = x i + 1 − x , a + b = h = Δ x , 则 6 h f i ( x ) = 6 ( a y i + 1 + b y i ) + a ( a 2 − h 2 ) M i + 1 + b ( b 2 − h 2 ) M i M i 为节点的二阶导数, 对应于力学上的弯矩, 满足下面的方程 M i + 4 M i + 1 + M i + 2 = d i + 1 = 6 h 2 ( y i + 2 − 2 y i + 1 + y i ) , i = 1 , 2... n − 2 要求的 M i 个数为 n , 而对应的方程数目为 n − 2 , 故还需两个边界条件才能唯一确定, 边界条件可取为两端点的导数值或是二阶导数值. 常用的自然边界条件指 M 1 = M n = 0 . 加上边界条件后便可求得 M 2 , M 3 , . . . M n − 1 M i 满足的方程为 0 M 2 M 3 ⋮ M n − 3 M n − 2 + + + ⋮ + + 4 M 2 4 M 3 4 M 4 ⋮ 4 M n − 2 4 M n − 1 + + + ⋮ + + M 3 M 4 M 5 ⋮ M n − 1 0 = = = ⋮ = = d 2 d 3 d 4 d n − 2 d n − 1 写为矩阵形式 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 4 1 ⋮ … … 1 4 ⋮ 1 0 0 1 ⋮ 4 1 … … ⋮ 1 4 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ M 2 M 3 ⋮ M n − 2 M n − 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ d 2 d 3 ⋮ d n − 2 d n − 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 可见, 此方程为三对角方程组, 对角占优, 存在唯一解, 可利用所谓的追赶法求解. 中文常称的追赶法, 是 Thomas方法 的形象翻译. 大致求解过程分为两步: 追: 利用消元法将原方程化为二对角方程, 向前递推, 使系数矩阵主对角线变为1. 由第一个方程, 得到 M 2 和 M 3 的关系, 将其带入第二个方程, 消去主对角线下方系数. 以此进行, 最终追到 M n − 1 , 得到其解. 变换后, 其方程为 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1 0 ⋮ … … A 2 1 ⋮ 0 0 0 A 3 ⋮ 1 0 … … ⋮ A n − 2 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ M 2 M 3 ⋮ M n − 2 M n − 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ D 2 D 3 ⋮ D n − 2 D n − 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 赶: M n − 1 已经追得, 然后由此倒推, 得到其他 M i 值. 对上面的方程, 由于系数是固定的, Thomas方法的递推式为 A 2 D 2 M n − 1 = 1 4 , = A 2 d 2 , = D n − 1 , A i D i M i = 1 4 − A i − 1 = A i ( d i − D i − 1 ) = D i − A i M i + 1 对 A i , 可求得其通式 A i = 2 − 3 √ t i + 1 t i − 1 , t = ( 2 + 3 √ ) 2 对 D i , 向后递推至 D 2 , 一般项可写为 D j = ∑ k = 2 j ( − 1 ) j − k P j k d k 对 M i , 向前递推至 M n − 1 , 一般项可写为 M i = ∑ j = i n − 1 ( − 1 ) j − i P j − 1 i D j 综合上面两个结果, 得到 M i P n m = ∑ j = i n − 1 ∑ k = 2 j ( − 1 ) 2 j − i − k P j − 1 i P j k d k = ∑ j = i n − 1 ∑ k = 2 j ( − 1 ) i + k P j − 1 i P j k d k , i = 2 , 3 , . . . n − 1 = ∏ l = m n A l = A m A m + 1 … A n − 1 A n , n m 根据插值公式 d k 6 h f i ( x ) = 6 h 2 ( y k + 1 − 2 y k + y k − 1 ) = 6 ( α y i + 1 + β y i ) + α ( α 2 − h 2 ) M i + 1 + β ( β 2 − h 2 ) M i 以划分间距 h 为单位, 约化上述公式 f i ( x ) α = α y i + 1 + β y i + α ( α 2 − 1 ) μ i + 1 + β ( β 2 − 1 ) μ i = a h , β = b h = 1 − α , μ i = M i 6 / h 2 由此可见, 虽然三次样条函数仍可写为节点值的线性形式, 但其系数十分复杂. 上面公式看起来清楚, 但是实际计算时需要计算 M i 和 M i + 1 , 两个三重循环, 整体计算量为 O ( N 3 ) . 利用 A i 的近似关系和 M i 的递推关系可以将公式的计算量减少一些. f i ( x ) = α y i + 1 + β y i + α ( α 2 − 1 ) μ i + 1 + β ( β 2 − 1 ) μ i = α y i + 1 + β y i + β ( β 2 − 1 ) D i + μ i + 1 对不等距划分, 上面的递推公式太过复杂, 很难写出一般项了. 样条函数可写为 f i ( x ) h i = a y i + 1 + b y i h i + a ( a 2 − h 2 ) M i + 1 + b ( b 2 − h 2 ) M i 6 h i = x i + 1 − x i , a = x − x i , b = h i − a M i 满足方程 h i M i + 2 ( h i + h i + 1 ) M i + 1 + h i + 1 M i + 2 = 6 ( y i + 2 − y i + 1 h i + 1 − y i + 1 − y i h i ) , i = 1 , 2 , … , n − 2 代码及测试测试结果 awk 实现的代码如下 awk ' BEGIN { Ndat = 0 } NF = = 3 { Ndat + + ; X = $ 2 ; Y = $ 3 } END { for ( i = 1 ; i Ndat ; i + + ) dX = X - X for ( i = 2 ; i Ndat ; i + + ) d = 6 * ( ( Y - Y ) / dX - ( Y - Y ) / dX ) A = 0 ; A = 0 D = 0 ; D = 0 for ( i = 2 ; i = Ndat - 1 ; i + + ) { ai = dX bi = 2 * ( dX + dX ) ci = dX A = ci / ( bi - ai * A ) D = ( d - ai * D ) * A / ci } M = 0 for ( i = Ndat - 1 ; i 1 ; i - - ) M = D - A * M # for ( i = 1 ; i = Ndat ; i + + ) print i , dX , d , A , D , M h = X - X for ( x = X ; x X ; x + = . 1 * h ) print x , SP 2 ( x , 0 , Ndat , X , Y , dX ) , SP 3 ( x , 0 , Ndat , X , Y , dX , M ) } function SP 2 ( x , i , Ndat , X , Y , dX , j , a , ki ) { if ( i = = 0 ) { for ( i = 1 ; i Ndat ; i + + ) if ( X x ) break } ki = 0 for ( j = i ; j = Ndat - 1 ; j + + ) ki + = 2 * ( - 1 ) ^ ( j - i ) * ( Y - Y ) / dX a = ( x - X ) / dX return a^ 2 * Y + ( 1 - a^ 2 ) * Y + ( 1 - a ) * ki * ( x - X ) } # function SP 3 ( x , i , Ndat , X , Y , dX , M , a , b , h ) { if ( i = = 0 ) { for ( i = 1 ; i Ndat ; i + + ) if ( X x ) break } h = dX a = x - X ; b = h - a return ( a * Y + b * Y ) / h + ( a * ( a^ 2 - h^ 2 ) * M + b * ( b^ 2 - h^ 2 ) * M ) / ( 6 * h ) } # function Psp 3 ( Ndat , i , j , k , a ) { a = ( 2 + sqrt ( 3 ) ) ^ 2 for ( i = 2 ; i Ndat ; i + + ) { for ( j = i ; j Ndat ; j + + ) { P = 1 for ( k = i ; k = j ; k + + ) P * = 2 - sqrt ( 3 ) * ( a^k + 1 ) / ( a^k - 1 ) } } } # function uniSP 3 ( x , i , Ndat , X , Y , j , k , a , b , h ) { if ( i = = 0 ) { for ( i = 1 ; i Ndat ; i + + ) if ( X x ) break } h = X - X a = ( x - X ) / h ; b = 1 - a for ( j = 1 ; j = Ndat ; j + + ) Coef = 0 Coef + = b Coef + = a Rsec = b * ( b^ 2 - 1 ) for ( k = 2 ; k = i ; k + + ) { Rtmp = Rsec * ( - 1 ) ^ ( i - k ) * P Coef + = Rtmp Coef + = - 2 * Rtmp Coef + = Rtmp } Rsec = a * ( a^ 2 - 1 ) - b * ( b^ 2 - 1 ) * P for ( j = i + 1 ; j = Ndat - 1 ; j + + ) { Rij = Rsec * ( - 1 ) ^ ( i + 1 ) if ( j ! = i + 1 ) Rij * = P for ( k = 2 ; k = j ; k + + ) { Rtmp = Rij * ( - 1 ) ^k * P Coef + = Rtmp Coef + = - 2 * Rtmp Coef + = Rtmp } } F = 0 for ( j = 1 ; j = Ndat ; j + + ) F + = Coef * Y return F # ( a * Y + b * Y + a * ( a^ 2 - 1 ) * Mi ( i + 1 , Ndat ) + b * ( b^ 2 - 1 ) * Mi ( i , Ndat ) ) } # function Mi ( i , Ndat , j , k , Rtmp , ret ) { ret = 0 for ( j = i ; j = Ndat - 1 ; j + + ) { Rtmp = ( - 1 ) ^i if ( j ! = i ) Rtmp * = P for ( k = 2 ; k = j ; k + + ) ret + = Rtmp * ( - 1 ) ^k * P * d } return ret } ' ABS ABS . sp # ' CUB CUB.sp #' LJ LJ.sp 利用Matlab的 csape 函数 可以进行三次样条函数的插值, 示例代码如下 format long x = 1 : 0 . 02 : 15 ; y = - 1 E 3 * ( - 12 . / x . ^ 13 + 6 . / x . ^ 7 ) ; pp = csape ( x , y , 'second' , ) ; pp = csape ( x , y ) ; xsp = 1 : 0 . 01 : 1 . 5 ; ysp = ppval ( pp , xsp ) ; yst = - 1 E 3 * ( - 12 . / xsp . ^ 13 + 6 . / xsp . ^ 7 ) ; plot ( x,y, 'o',xsp,ysp,'-',xsp,yst,'-' ) axis ( ) FID = fopen ( 'LJ.mat', 'w' ) ; Ndat = length ( xsp ) ; for i = 1 : Ndat fprintf ( FID , '%f %f\n' , xsp ( i ) , ysp ( i ) ) ; end 几个测试函数的结果 样条函数的积分 对已经拟合的样条函数进行数值积分时可利用可利用 Simpson方法 . Simpson方法具有三阶精度, 对不超过三次的多项式精确成立, 很适合于二次和三次样条函数.
以前总结的复杂网络和人类动力学研究中常见的分布律及拟合方法，挂在博客上方便取用，也欢迎各位同行指教！ 目录： 基本术语英汉对照 常见的分布律 l 正态 / 高斯分布 Normal distribution / Gaussian distribution l 对数正态分布 Log-normal distribution l 指数 / 负指数分布 Exponential distribution / Negative exponential distribution l 泊松分布 Poisson distribution l 幂律分布 Power law distribution l 指数截断的幂律分布 Power law with exponential cutoff l 截断幂律 Truncated power law l 广延指数分布 Stretched exponential distribution l 漂移幂律 Shifted power law 数据拟合与参数估计（7步） 重要参考文献 常见分布律及数据拟合总结.docx
gnuplot 除了绘图功能之外，最简单实用的功能就是拟合了。gnuplot 可以进行单变量甚至多变量的线性和非线性拟合。虽然可能不像专门的数学软件那么强大，但是足以对付日常需要了。我们拿上一篇文章里的数据来举例子。 首先，要定义一个待拟合的函数： gnuplot f(x)=50*(1+erf(a*(x-b))) 这里使用了误差函数 erf(x) ，有两个待定的参数： a, b 。下面我们生成一个文件“ fit.par ”，里面包含的是参数 a 和 b 的初值： a = 1 b = 12 初值的选择要尽可能贴近结果，否则可能导致误差甚至无法收敛。下面我们进行拟合： gnuplot fit f(x) 'probability.dat' using 1:2 via 'fit.par' gnuplot 里面关于拟合的命令是 fit ，后面的自变量取值范围不是必需的。 f(x) 函数已经在上面定义过了，数据文件“ probability.dat ”也已经在上一篇博文中交代过了。 via 后面跟的是参数变量列表文件。执行 fit 命令之后，gnuplot 会输出一堆结果。我们忽略那些中间运算，只把最后结果贴在下面： After 5 iterations the fit converged. final sum of squares of residuals : 41.9399 rel. change during last iteration : -4.27973e-07 degrees of freedom (FIT_NDF) : 8 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 2.28965 variance of residuals (reduced chisquare) = WSSR/ndf : 5.24249 Final set of parameters Asymptotic Standard Error ======================= ========================== a = 1.15661 +/- 0.06331 (5.474%) b = 11.9027 +/- 0.02383 (0.2002%) correlation matrix of the fit parameters: a b a 1.000 b 0.014 2.000 这段文字说明，经过 5 次迭代，gnuplot 得到了收敛的结果。中间部分是参数 a 和 b 的最终取值以及渐近标准差（asymptotic standard error）。渐近标准差的计算是基于线性拟合的，对于非线性拟合，渐近标准差一般都比真的标准差小，所以这个数字只能用于定性分析。而最后给出的相关矩阵（correlation matrix）可以帮助我们确认渐近标准差的可靠度，非对角元素绝对值越小，渐近标准差越接近真实标准差。 好了，现在我们可以把数据和拟合曲线画在同一张图上了： gnuplot set xrange gnuplot set yrange gnuplot unset key gnuplot set xlabel "Laser Pulse Energy (μJ)" gnuplot set ylabel "Bubble Formation Probability (%)" gnuplot plot "probability.dat" using 1:2:3:4 with xerrorbars, f(x) lw 2 lc rgb "orange"
标签: 拟合
