《天文数据分析课件》由会员分享，可在线阅读，更多相关《天文数据分析课件（50页珍藏版）》请在金锄头文库上搜索。
1、天文数据分析,国家天文台 赵永恒 2015年4月,数据分析方法,累积概率分布,任意分布的随机数,常用概率分布,数据插值,内插和外推,插值(scipy.interpolate),线性1d插值 (interp1d) 1d样条插值 (interpolate.splXXX) 2d样条插值 (bisplrep) RBF（radial basis function）插值/平滑,import numpy as np from scipy.interpolate import Rbf import matplotlib.pyplot as plt from matplotlib import cm # 2-d
2、 tests - setup scattered data x = np.random.rand(100)*4.0-2.0 y = np.random.rand(100)*4.0-2.0 z = x*np.exp(-x*2-y*2) ti = np.linspace(-2.0, 2.0, 100) XI, YI = np.meshgrid(ti, ti) # use RBF rbf = Rbf(x, y, z, epsilon=2) ZI = rbf(XI, YI) # plot the result n = plt.Normalize(-2., 2.) plt.subplot(1, 1, 1
3、) plt.pcolor(XI, YI, ZI, cmap=cm.jet) plt.scatter(x, y, 100, z, cmap=cm.jet) plt.title(RBF interpolation - multiquadrics) plt.xlim(-2, 2) plt.ylim(-2, 2) plt.colorbar() plt.show(),数据拟合,拟合模型,O = T + e,最小二乘法,import numpy as np import scipy.linalg as splin import matplotlib.pyplot as plt c1,c2= 5.0,2.0
4、 i = np.r_1:11 xi = 0.1*i yi = c1*np.exp(-xi)+c2*xi zi = yi + 0.05*np.max(yi)*np.random.randn(len(yi) A = np.c_np.exp(-xi):,np.newaxis,xi:,np.newaxis c,resid,rank,sigma = splin.lstsq(A,zi) print c,resid,rank,sigma xi2 = np.r_0.1:1.0:100j yi2 = c0*np.exp(-xi2) + c1*xi2 plt.plot(xi,zi,x,xi2,yi2) plt.a
5、xis(0,1.1,3.0,5.5) plt.xlabel($x_i$) plt.title(Data fitting with linalg.lstsq) plt.show(),非线性拟合,拟合系数有非线性函数关系,from numpy import x = arange(0,6e-2,6e-2/30) A,k,theta = 10, 1.0/3e-2, pi/6 y_true = A*sin(2*pi*k*x+theta) y_meas = y_true + 2*random.randn(len(x) def residuals(p, y, x): A,k,theta = p err =
6、y-A*sin(2*pi*k*x+theta) return err def peval(x, p): return p0*sin(2*pi*p1*x+p2) p0 = 8, 1/2.3e-2, pi/3 print array(p0) # 8. 43.4783 1.0472 from scipy.optimize import leastsq plsq = leastsq(residuals, p0, args=(y_meas, x) print plsq0 # 10.9437 33.3605 0.5834 print array(A, k, theta) # 10. 33.3333 0.5
7、236 import matplotlib.pyplot as plt plt.plot(x,peval(x,plsq0),x,y_meas,o,x,y_true) plt.title(Least-squares fit to noisy data) plt.legend(Fit, Noisy, True) plt.show(),Maximum Likelihood Estimation (MLE),最大似然估计,一致收敛 渐近于正态分布：有极值,等精度高斯分布,非等精度高斯分布,上下限分布,求极值scipy.optimize,A collection of general-purpose o
8、ptimization routines. fmin - Nelder-Mead Simplex algorithm (uses only function calls) fmin_powell - Powells (modified) level set method (uses only function calls) fmin_cg - Non-linear (Polak-Ribiere) conjugate gradient algorithm (can use function and gradient). fmin_bfgs - Quasi-Newton method (Broyd
9、on-Fletcher-Goldfarb-Shanno); (can use function and gradient) fmin_ncg - Line-search Newton Conjugate Gradient (can use function, gradient and Hessian). leastsq - Minimize the sum of squares of M equations in N unknowns given a starting estimate. Constrained Optimizers (multivariate) fmin_l_bfgs_b -
10、 Zhu, Byrd, and Nocedals L-BFGS-B constrained optimizer (if you use this please quote their papers - see help) fmin_tnc - Truncated Newton Code originally written by Stephen Nash and adapted to C by Jean-Sebastien Roy. fmin_cobyla - Constrained Optimization BY Linear Approximation,EM算法,Expect-Maximu
11、m算法 EM算法的目标是找出有隐性变量的概率模型的最大可能性解 分为2个步骤，E-step和M-step E-step根据最初假设的模型参数值或者上一步的模型参数计算出隐性变量的后验概率，其实就是隐性变量的期望 M-step根据这个E-step的后验概率重新计算出模型参数，然后再重复这两个步骤，直至目标函数收敛。,拟合度,好模型,精度估计,Bootstrap方法 “杰克刀”方法,Bootstrap方法,用已知数据建立分布 可有N!组新数据,SN1987A中微子事件,9个中微子,“杰克刀”方法,每次丢掉1或多个数据,KolmogorovSmirnov (K-S) 检验,高斯分布？,Anderso
12、nDarling检验 ShapiroWilk检验,Bayes统计推断,高斯分布,问题：,问题：,直方图,Markov chain Monte Carlo (MCMC),https:/ 5-1,LombScargle 周期图,不等间隔,相位图,注意：倍频、差频,相关分析,小波变换,课后作业II,1、一组有确定均值和误差的测量值中，即有N个值xi，又有M个上下限值xmini,xmaxi，请写出其似然函数ln L。 2、随机写出20个整数（可重复），组成xi。以xi的分布产生出1000个随机数，构成yi。并用KS检验说明xi和yi的分布是否一致。 3、设计具有一个或数个频率（周期）的正弦曲线，并加上误差。用Fourier谱分析确定其频率（周期）。,
