loess局部平滑

代码

#-*-coding:gbk-*-
#########################################################################
#   Copyright (C) 2013 All rights reserved.
# 
#   文件名称：loess.py
#   创 建 者：刘禹 finallyly 
#   创建日期：2013年01月27日
#   描    述：
#
#   备    注：
#
#########################################################################
#!/usr/bin/python
# please add your code here!

import re;
import sys;
import time;
import matplotlib
matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
from pylab import *
# x location; h:bandwith;xp,yp:data points(vectors)
def loess(x,h,xp,yp):
    w=exp(-0.5*((x-xp)/h)**2)/sqrt(2*pi*h**2);
    b=sum(w*xp)*sum(w*yp)-sum(w)*sum(w*xp*yp);
    b/=sum(w*xp)**2-sum(w)*sum(w*xp**2);
    a=(sum(w*yp)-b*sum(w*xp))/sum(w);
    return a+b*x;
def PrintUsage():
    sys.stderr.write("program [IN]\n");
    exit(1);
if(__name__=="__main__"):

    if(len(sys.argv)!=2):
        PrintUsage();
    starttime=time.clock();
    d=loadtxt(str(sys.argv[1]));
    s1=[];
    s2=[];  
    for k in d[:,0]:
        s1.append(loess(k,5,d[:,0],d[:,1]));
        s2.append(loess(k,100,d[:,0],d[:,1]));
    xlabel("Day in Year");
    ylabel("Draf Number");
    gca().set_aspect('equal');
    plot(d[:,0],d[:,1],'o',color="white",markersize=7,linewidth=3);
    plot(d[:,0],s1,'k-',d[:,0],s2,'k--');
    q=4;
    axis([1-q,366+q,1-q,366+q]);
    savefig("draftlottery.png");
    endtime=time.clock();
    interval=endtime-starttime;
    sys.stderr.write("%s has finished congratulations!\n"%str(sys.argv[0]));
    sys.stderr.write("time elapse:%f\n"%interval);

 数据地址 http://www.amstat.org/publications/jse/v5n2/datasets.starr.html#rosenbaum1，下载draft70yr.dat.txt，去掉第一列空白

运行../python2.7 loess.py  draft70yr.dat.txt

结果：