天文数据分析课件

天文数据分析国家天文台赵永恒2015年4月天文数据分析国家天文台赵永恒1数据分析方法数据分析方法2累积概率分布累积概率分布3任意分布的随机数任意分布的随机数4常用概率分布常用概率分布5数据插值数据插值6内插和外推内插和外推7插值(erpolate)线性1d插值(interp1d)1d样条插值(interpolate.splXXX)2d样条插值(bisplrep)RBF（radialbasisfunction）插值/平滑插值(erpolate)线性1d插值(i8importnumpyasnpfromerpolateimportRbfimportmatplotlib.pyplotaspltfrommatplotlibimportcm#2-dtests-setupscattereddatax=np.random.rand(100)*4.0-2.0y=np.random.rand(100)*4.0-2.0z=x*np.exp(-x**2-y**2)ti=np.linspace(-2.0,2.0,100)XI,YI=np.meshgrid(ti,ti)#useRBFrbf=Rbf(x,y,z,epsilon=2)ZI=rbf(XI,YI)#plottheresultn=plt.Normalize(-2.,2.)plt.subplot(1,1,1)plt.pcolor(XI,YI,ZI,cmap=cm.jet)plt.scatter(x,y,100,z,cmap=cm.jet)plt.title('RBFinterpolation-multiquadrics')plt.xlim(-2,2)plt.ylim(-2,2)plt.colorbar()plt.show()importnumpyasnp9数据拟合数据拟合10拟合模型O=T+e拟合模型O=T+e11最小二乘法最小二乘法12importnumpyasnpimportscipy.linalgassplinimportmatplotlib.pyplotaspltc1,c2=5.0,2.0i=np.r_[1:11]xi=0.1*iyi=c1*np.exp(-xi)+c2*xizi=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)printc,resid,rank,sigmaxi2=np.r_[0.1:1.0:100j]yi2=c[0]*np.exp(-xi2)+c[1]*xi2plt.plot(xi,zi,'x',xi2,yi2)plt.axis([0,1.1,3.0,5.5])plt.xlabel('$x_i$')plt.title('Datafittingwithlinalg.lstsq')plt.show()importnumpyasnp13非线性拟合拟合系数有非线性函数关系非线性拟合拟合系数有非线性函数关系14fromnumpyimportx=arange(0,6e-2,6e-2/30)A,k,theta=10,1.0/3e-2,pi/6y_true=A*sin(2*pi*k*x+theta)y_meas=y_true+2*random.randn(len(x))defresiduals(p,y,x):A,k,theta=perr=y-A*sin(2*pi*k*x+theta)returnerrdefpeval(x,p):returnp[0]*sin(2*pi*p[1]*x+p[2])p0=[8,1/2.3e-2,pi/3]printarray(p0)#[8.43.47831.0472]fromscipy.optimizeimportleastsqplsq=leastsq(residuals,p0,args=(y_meas,x))printplsq[0]#[10.943733.36050.5834]printarray([A,k,theta])#[10.33.33330.5236]importmatplotlib.pyplotaspltplt.plot(x,peval(x,plsq[0]),x,y_meas,'o',x,y_true)plt.title('Least-squaresfittonoisydata')plt.legend(['Fit','Noisy','True'])plt.show()fromnumpyimport15MaximumLikelihoodEstimation(MLE)MaximumLikelihoodEstimation16最大似然估计一致收敛渐近于正态分布：有极值最大似然估计一致收敛17等精度高斯分布等精度高斯分布18非等精度高斯分布非等精度高斯分布19上下限分布上下限分布20求极值scipy.optimizeAcollectionofgeneral-purposeoptimizationroutines.fmin--Nelder-MeadSimplexalgorithm(usesonlyfunctioncalls)fmin_powell--Powell's(modified)levelsetmethod(usesonlyfunctioncalls)fmin_cg--Non-linear(Polak-Ribiere)conjugategradientalgorithm(canusefunctionandgradient).fmin_bfgs--Quasi-Newtonmethod(Broydon-Fletcher-Goldfarb-Shanno);(canusefunctionandgradient)fmin_ncg--Line-searchNewtonConjugateGradient(canusefunction,gradientandHessian).leastsq--MinimizethesumofsquaresofMequationsinNunknownsgivenastartingestimate.ConstrainedOptimizers(multivariate)fmin_l_bfgs_b--Zhu,Byrd,andNocedal'sL-BFGS-Bconstrainedoptimizer(ifyouusethispleasequotetheirpapers--seehelp)fmin_tnc--TruncatedNewtonCodeoriginallywrittenbyStephenNashandadaptedtoCbyJean-SebastienRoy.fmin_cobyla--ConstrainedOptimizationBYLinearApproximation求极值scipy.optimizeAcollection21EM算法Expect-Maximum算法EM算法的目标是找出有隐性变量的概率模型的最大可能性解分为2个步骤，E-step和M-stepE-step根据最初假设的模型参数值或者上一步的模型参数计算出隐性变量的后验概率，其实就是隐性变量的期望M-step根据这个E-step的后验概率重新计算出模型参数，然后再重复这两个步骤，直至目标函数收敛。EM算法Expect-Maximum算法22拟合度拟合度23好模型好模型24精度估计Bootstrap方法“杰克刀”方法精度估计Bootstrap方法25Bootstrap方法用已知数据建立分布可有N!组新数据Bootstrap方法用已知数据建立分布26SN1987A中微子事件9个中微子SN1987A中微子事件9个中微子27“杰克刀”方法每次丢掉1或多个数据“杰克刀”方法每次丢掉1或多个数据28Kolmogorov–Smirnov(K-S)检验Kolmogorov–Smirnov(K-S)检验29高斯分布？Anderson–Darling检验Shapiro–Wilk检验

高斯分布？Anderson–Darling检验30Bayes统计推断Bayes统计推断31高斯分布问题：高斯分布问题：32问题：问题：33天文数据分析课件34天文数据分析课件35直方图直方图36MarkovchainMonteCarlo(MCMC)/pymc-devs/pymc

MarkovchainMonteCarlo(MCMC37时间序列分析时间序列分析38Fourier变换scipy.fftpackFourier变换scipy.fftpack39卷积定理卷积定理40离散Fourier变换离散Fourier变换41窗效应窗效应42GX5-1GX5-143Lomb–Scargle周期图不等间隔Lomb–Scargle周期图不等间隔44相位图注意：倍频、差频…相位图注意：倍频、差频…45相关分析相关分析46小波变换小波变换47天文数据分析课件48天文数据分析课件49课后作业II1、一组有确定均值和误差σ的测量值中，即有N个值{xi}，又有M个上下限值{xmini,xmaxi}，请写出其似然函数lnL。2、随机写出20个整数（可重复），组成{xi}。以{xi}的分布产生出1000个随机数，构成{yi}。并用KS检验说明{xi}和{yi}的分布是否一致。3、设计具有一个或数个频率（周期）的正弦曲线，并加上误差。用Fourier谱分析确定其频率（周期）。课后作业II1、一组有确定均值和误差σ的测量值中，即有N个50天文数据分析国家天文台赵永恒2015年4月天文数据分析国家天文台赵永恒51数据分析方法数据分析方法52累积概率分布累积概率分布53任意分布的随机数任意分布的随机数54常用概率分布常用概率分布55数据插值数据插值56内插和外推内插和外推57插值(erpolate)线性1d插值(interp1d)1d样条插值(interpolate.splXXX)2d样条插值(bisplrep)RBF（radialbasisfunction）插值/平滑插值(erpolate)线性1d插值(i58importnumpyasnpfromerpolateimportRbfimportmatplotlib.pyplotaspltfrommatplotlibimportcm#2-dtests-setupscattereddatax=np.random.rand(100)*4.0-2.0y=np.random.rand(100)*4.0-2.0z=x*np.exp(-x**2-y**2)ti=np.linspace(-2.0,2.0,100)XI,YI=np.meshgrid(ti,ti)#useRBFrbf=Rbf(x,y,z,epsilon=2)ZI=rbf(XI,YI)#plottheresultn=plt.Normalize(-2.,2.)plt.subplot(1,1,1)plt.pcolor(XI,YI,ZI,cmap=cm.jet)plt.scatter(x,y,100,z,cmap=cm.jet)plt.title('RBFinterpolation-multiquadrics')plt.xlim(-2,2)plt.ylim(-2,2)plt.colorbar()plt.show()importnumpyasnp59数据拟合数据拟合60拟合模型O=T+e拟合模型O=T+e61最小二乘法最小二乘法62importnumpyasnpimportscipy.linalgassplinimportmatplotlib.pyplotaspltc1,c2=5.0,2.0i=np.r_[1:11]xi=0.1*iyi=c1*np.exp(-xi)+c2*xizi=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)printc,resid,rank,sigmaxi2=np.r_[0.1:1.0:100j]yi2=c[0]*np.exp(-xi2)+c[1]*xi2plt.plot(xi,zi,'x',xi2,yi2)plt.axis([0,1.1,3.0,5.5])plt.xlabel('$x_i$')plt.title('Datafittingwithlinalg.lstsq')plt.show()importnumpyasnp63非线性拟合拟合系数有非线性函数关系非线性拟合拟合系数有非线性函数关系64fromnumpyimportx=arange(0,6e-2,6e-2/30)A,k,theta=10,1.0/3e-2,pi/6y_true=A*sin(2*pi*k*x+theta)y_meas=y_true+2*random.randn(len(x))defresiduals(p,y,x):A,k,theta=perr=y-A*sin(2*pi*k*x+theta)returnerrdefpeval(x,p):returnp[0]*sin(2*pi*p[1]*x+p[2])p0=[8,1/2.3e-2,pi/3]printarray(p0)#[8.43.47831.0472]fromscipy.optimizeimportleastsqplsq=leastsq(residuals,p0,args=(y_meas,x))printplsq[0]#[10.943733.36050.5834]printarray([A,k,theta])#[10.33.33330.5236]importmatplotlib.pyplotaspltplt.plot(x,peval(x,plsq[0]),x,y_meas,'o',x,y_true)plt.title('Least-squaresfittonoisydata')plt.legend(['Fit','Noisy','True'])plt.show()fromnumpyimport65MaximumLikelihoodEstimation(MLE)MaximumLikelihoodEstimation66最大似然估计一致收敛渐近于正态分布：有极值最大似然估计一致收敛67等精度高斯分布等精度高斯分布68非等精度高斯分布非等精度高斯分布69上下限分布上下限分布70求极值scipy.optimizeAcollectionofgeneral-purposeoptimizationroutines.fmin--Nelder-MeadSimplexalgorithm(usesonlyfunctioncalls)fmin_powell--Powell's(modified)levelsetmethod(usesonlyfunctioncalls)fmin_cg--Non-linear(Polak-Ribiere)conjugategradientalgorithm(canusefunctionandgradient).fmin_bfgs--Quasi-Newtonmethod(Broydon-Fletcher-Goldfarb-Shanno);(canusefunctionandgradient)fmin_ncg--Line-searchNewtonConjugateGradient(canusefunction,gradientandHessian).leastsq--MinimizethesumofsquaresofMequationsinNunknownsgivenastartingestimate.ConstrainedOptimizers(multivariate)fmin_l_bfgs_b--Zhu,Byrd,andNocedal'sL-BFGS-Bconstrainedoptimizer(ifyouusethispleasequotetheirpapers--seehelp)fmin_tnc--TruncatedNewtonCodeoriginallywrittenbyStephenNashandadaptedtoCbyJean-SebastienRoy.fmin_cobyla--ConstrainedOptimizationBYLinearApproximation求极值scipy.optimizeAcollection71EM算法Expect-Maximum算法EM算法的目标是找出有隐性变量的概率模型的最大可能性解分为2个步骤，E-step和M-stepE-step根据最初假设的模型参数值或者上一步的模型参数计算出隐性变量的后验概率，其实就是隐性变量的期望M-step根据这个E-step的后验概率重新计算出模型参数，然后再重复这两个步骤，直至目标函数收敛。EM算法Expect-Maximum算法72拟合度拟合度73好模型好模型74精度估计Bootstrap方法“杰克刀”方法精度估计Bootstrap方法75