Lattice复习笔记I

Lattice这个词貌似是一个GRE词汇,原意是指格子框架之类的东西,而在R中则是指Lattice包。Lattice是一个多元数据绘图的包,同时为了完成绘图数据格式的要求,也内置了一些数据整理的函数。我最初接触这个包的时候,还需要从install.package中下载,但是后来就变成和MASS、Grid一样的自带包了,由此可见其“江湖地位”。但是令我比较不解的是,在我上研到毕业工作之四年中,R在国内的应用可以说是飞速发展,毕业两年基本没有再接触过R编程的我,前些日子无意间搜了一下Lattice包的相关文献,却发现仍然没有中文相关的文章来系统地介绍这个包,只有几篇简单的轮廓性描写的文章,甚至没有一篇“自主创新”的应用类文章。于是我萌生了写一个Lattice包相关的系列文章的想法,工作养成的习惯,有了大致思路后就开始上手,發現 IQ Option 但是细节与系统性方面可能有些欠缺,还望读者见谅,主要的参考文献,还是包的作者Deepayan Sarkar的那本《Lattice: Multivariate Data Visualization with R》。

The ‘lattice’ add-on package is an implementation of Trellis graphics for R. It is a powerful and elegant high-level data visualization system with an emphasis on multivariate data. It is designed to meet most typical graphics needs with minimal tuning, but can also be easily extended to handle most nonstandard requirements.

还是从最基本的介绍开始吧,在R中输入?lattice后出现三段话,第一段是包的基本描述,这段话里有几个关键词,一是Trellis graphics,也就是格子作图;二是Multivariate data,即多元数据作图;三是nonstandard requirements。看起来很抽象,画出来就容易理解了。

lattice_1

所谓格子作图,可以直观地理解为“格子集合”,每一个“格子”都包含一张独立的统计图形。上面的样图可以视作四个格子按2乘2矩阵格式排列而成的结果,只不过(1,2)位置是空白。多元数据指的是多变量的数据,比如上面作图所用的数据是R里自带的iris数据集,这个数据集经常被当作判别分析的例子,使用IQ Option 它包括五个变量,如下图:

image

四个数值型变量和一个分类型变量,实际上,如果使用lattice包中的equal.count函数的话,可以将连续型的数据变量分割成不同的数据区间(Lattice中称为shingles),从而将其转化为分类变量(每个数据区间即为一类)。比如我们可以把Sepal.Length变量分为大于5和小于5两类,然后就可以将其当作分类型变量来使用了。而nonstandard requirements指的是lattice包中提供了一些基础的作图函数,而这些函数又有多种的选项,这些选项可以满足用户定制化的需求。比如你想将上图中的格子全部放在一行中,或者想将直方图的柱子换个颜色等等,就可以通过改变对应的参数来完成。

High-level ‘lattice’ functions like ‘xyplot’ are different from traditional R graphics functions in that they do not perform any plotting themselves.  Instead, they return an object, of class ‘”trellis”‘, which has to be then ‘print’-ed or ‘plot’-ted to create the actual plot.  Due to R’s automatic printing rule, it is usually not necessary to explicitly carry out the second step, and ‘lattice’ functions appear to behave like their traditional counterparts.  However, the automatic plotting is suppressed when the high-level functions are called inside another function (most often ‘source’) or in other contexts where automatic printing is suppressed (e.g., ‘for’ or ‘while’ loops).  In such situations, an explicit call to ‘print’ or ‘plot’ is required.

第二段话是详细介绍,后面我们会作说明,现在直接看第三段话。这段话主要说明了lattice包与R传统作图的区别,也就是lattice的作图过程。实际上,所有lattice包内的作图函数,本身是不产生任何图形的,它们只会返回一个叫“Trellis”类型的对象,只有当这个对象被print或者plot方法调用时,才会产生出可视化的图形。因此lattice的作图分为两步,第一步是返回Trellis类型的对象,第二步是被print、plot过程调用,生成图形。只不过一般环境下第二步被默认执行,因此整个过程看起来像是一步完成的。注意下面的代码,

image

当x = hist()一行代码完成后,回车即显示图形,但是当y = histogram()一行完成时,没有任何图形产生,只有执行plot函数图形才会出现。类似的环境还有被其他通用函数调用时,或者在for和while的循环中时,lattice都是不会直接生成图形的。使用class函数分别查看x和y的类型,x为histogram,而y则是trellis类。如果想要对y进行某些参数上的调整,可以使用update函数,例如update(y, col = “blue”)可以将直方图调整为蓝色填充。而使用update(y, col = “blue”, layout = c(1, 3))则不仅可以调整图形的颜色,还可以调整其布局,如下图。

lattice_2

本小节暂时结束,下一小节我将先介绍lattice中作图函数的基本结构,然后开始总结lattice中的各种参数的调整方法,工作量不小,希望能坚持写完这个系列的文章。

美股开户过程记录

江堂兄的博客里有一篇专门介绍如何申请美国本土工作机会的文章,我记得里面他提到过申请美国工作最大的问题可能是,你根本就没想过去美国工作。言外之意,只要你有这个意识再加上一定的执行力和决心,去美国工作也并不是一件天方夜谭的事。其实申请美股账户也一样,最大困难是你根本没想过开美股的账户。而二者不同的地方在于,去美国工作要受到自身很多条件的限制,并不是一件100%能做成的事情,而开美股账户,只要有执行力,就会是一个必然能成的事情。

开户的第一件事是选券商,美股的券商很多,但是国内的人,尤其新手的话一般都在Sogotrade、Firstrade和Scottrade三者之间选择。单从交易费用上来说的话,Sogotrade是最低的,每笔交易(单向)大约是3美金,而Firstrade和Scottrade分别是6.95和7美金,Sogotrade有很大的优势,期权价格也大体是这个情况。这里不做详述,网上有很多相关的对比文章,随便搜一下就可以了。我最后选的是Firstrade,主要原因有三,一是本身我的交易频率并不高,交易费的支出也相应较少,二是Firstrade可以购买美股市场上种类繁多的基金,而其他两家不行,三是因为Firstrade整个的开户体验非常好,当然这是开户之后才知道的。

image

登录Firstrade的官网,点注册,开始流程化的填写,按自身的条件如实填写即可,需要注意的有两点:

1,账户类型分两种,一种是现金(Cash)账户,一种是融资(Margin)账户,顾名思义前者是只能用自己的钱来购买股票,而后者则可以借券商的钱来购买(年利率约7.5%,月结)。现金账户没有最低资金要求限制,但是实行T+3,直观理解是当天买过的股票必须三天后才能卖出,否则账户会被判违规,90天内不能交易,而融资账户则有2000刀的最低资金要求,且5天内有三次T+0的机会。许多人认为美股完全实行T+0,其实是不对的。

2,个人的家庭住址直接将实际住址翻译成英文即可,当然“路”、“街道”什么的还是得写成Road,Street,Firstrade不要求住址证明文件(水电缴费单、信用卡账单什么的),而有的券商是要求的。

注册完成后,邮箱里会收到一封邮件,要求提交材料。

image

注意邮件里的开户文件四个字是超链接,点开后可以登录账户,下载那两份文件的PDF版,打印出来签好名字,签名用中文,日期按要求格式填。护照的复印件,个人资料页和签名在同一页上,扫描之即可。

将这三份文件扫描成PDF文件后,电邮至Newaccounts@firstrade.com,注意必须是扫描,拍照不行。大概过一到两个工作日后,邮箱里会收到一封邮件,提示开户成功。

image

之前我在网上看的教程都说需要邮寄原版纸质的W8Ben文件到美国,为此我还特地咨询了Firstrade,得到的答复是他们不需要任何纸质文件,不需要邮寄文件过去。

开户完成后就到了最后一步,向账户里注入资金,一般的方法是电汇过去,或者直接开招行的香港账户(在北京就可以完成,网上有教程)。之所以要绕这么个弯,是因为国家目前不允许个人向境外机构进行投资。当然别怕,最坏的结果也只是钱被退回来(不知道手续费是否会被扣)。

我直接用工行的网银完成的汇款,前提是需要有个U盾,步骤如下:

1,在“结售汇”一项下,选择“购汇”,汇钞标志选默认的“汇”即可(千万别改选“钞”),币种是美元,用途选“旅游”即可;

2,购汇完成后,在转账汇款中选“跨境汇款-向境外银行汇款”,按Firstrade里面提供的收款人信息一顿复制粘贴之后,信息填写完毕。这里要特别注意的是,在汇款人附言里,只写自己的Firstrde数字账户号和英文名字,其他一概不用填写,而汇款的目的选旅游。在复制收款人地址的时候,系统会提示什么街道不准之类的,不用理睬。

工行汇款的手续费是千分之0.8,电报费100元,有一定星级的话可以打个折扣。之后就是等待的过程,由于跨境汇款可能会经同某些中间银行,他们也会收取一定的手续费,因此如果你想恰好汇2000,5000等整数资金时,最好再多加上几十刀。我的汇款大约用了四天时间完成,在此期间可以每天网银查询一下汇款进度,当汇款跟踪显示“银行处理成功”的时候,资金就是快到账了。到账之后,Firstrade也会有提示的邮件。

image

至此就可以正式开始你的美股之旅了,看行情的软件我强力推荐腾讯的“自选股”,我个人觉得客户体验非常好,平时也可以多泡一泡雪球,学习一些投资的经验。

最后还要提示一句,股市有风险,入市需谨慎!

从WordPress博客到微信公众号

曾几何时,博客还是网络上非常流行的平台,大概是我上大学的2005年,新浪、网易、腾讯等各大门户网站都推出了自己的博客平台,一时间百花竞放,蔚为壮观。但是没过多久,校内网开始流行,之后等到我读研究生的时候,网络又成了微博的天下,等到毕业了,微信又成了最火的名词,玩似的就把飞信变成了历史名词,晓伟不禁由衷地感叹,不是我不明白,这世界变化太快,若想不被淘汰,必需紧跟时代。

自从四个月前写了两篇文章后,晓伟的博客又暂停了更新,这个习惯不好,又到了呼唤毅力,召唤执行力的时候。其实这四个月间晓伟并没有偷懒,也读了几本书,按道理来讲写几篇读书笔记还是不难的。但是问题出在我读的书上,大致是这么几本,章怡和的《往事并不如烟》、高华的《敏感词是如何升起的》以及麦克法夸尔的《敏感词的起源卷一》,从主流价值观的角度来看,这些书都是“反动书”中的战斗机,写这些书的读后感,会出现两大难点。一,书中牵涉的人物、事件太过于复杂,没有多年的研究功底,是绝不敢胡乱造次的;二,虽然晓伟博客属于玩票性质,点击量小,但是也害怕万一被举报而牵连了同IP上的其他网站。

前言扯得太长,本博客还是要传播积极的人生态度,因为人生本来不易,何必让残酷的历史真相平添压抑呢?言归正传,前不久微信更新了新版本,增强了自身的搜索功能,主界面的搜索框不仅能够搜本机的信息,而且能搜到网络中的公众号以及相关文章了。这个功能还是非常威武的,想想会不会有很多从前做网页推广的网站,逐渐地转向微信公众号的推广呢?这个搜索框背后的上亿用户,会带来多么恐怖的流量啊。受这点的启发,我突然想到,其实我这个写了四年的博客,文章大概有几十篇,除了一部分技术文章外,绝大多数都比较适合移动端的阅读,因此我为什么不建一个微信公众号呢?

说建就建,这就是传说中的执行力。上网一搜索之后,发现这个事其实不是很复杂,已经有很多现成的解决方案了,之前我还担心Wordpress文章向微信转移的问题,结果发现只要借助于一个插件,这二者就可以做到“无缝结合”。

1,登录微信公共平台,点击右上角的注册链接,按要求填写邮箱号及密码,随后邮箱里会收到一封确认邮件,点击邮件中的链接激活邮箱进入下一步;

2,选择是建立订阅号还是服务号,个人用户的话,选默认的订阅号即可,之后会有一个长长的信息登记页以及公众号信息页,按要求填写即可,需要注意的是公众号的名称并不是唯一的,是可以和现有公众号取相同名字的,比如你在公共号里搜王晓伟,就会看到两个相同名字的号,一个是我,另一个哥们估计我们同名;另一点需要注意的是,上传个人身份证照片时,这个照片一定要本人拿着身份证来拍照,照片中可以看到本人真实头像和身份证上的头像;

3,信息填写完毕提交后,系统提示会在7个工作日内审核完毕,根据我的情况,这个审核过程只有一到两个工作日,非常效率;

4,信息审核通过后就可以登录进自己的公共平台号了,在左边侧栏的公众号设置中完善一下相关的设置,然后在安全助手中开通相关的安全措施,比如手机扫描二维码登录可以防止公共号被盗用;

image

5,如果需要与Wordpress博客连接的话,可以在这里下载一个名为weixinpress的插件,在Wordpress的插件库中搜索不到这个插件,只能从上面的网址中下载该插件,然后手动安装,手动安装很容易,在插件页面中选择上传,然后选中下载的压缩文件即可;另外一个常见的的插件叫微信机器人,这个插件分为免费版和付费的高级版,免费版中没有自定义常见命令的功能,所以本文没有选用它;

image

6,插件安装好之后,进入设置页面,将最新文章、热门文章等常见功能的命令定义为1、2、3、4(初始值分别为New、Hot、Rand),因为考虑到英语单词毕竟还是太长,数字相对简洁明了许多;

image

7,进入微信公众平台的后端,进入设置中的开发者中心,启用服务器配置,在URL中填写你的博客地址,在Token中填写图6中weixinpress中红圈中的名称,要注意图6、7中红圈地方的的名称必须相同,否则微信将无法正常回应订阅者的命令;

image

8,保存设置后,在微信里搜索自己的公众号,然后关注,就可以看到欢迎界面了;

Welcome

回复命令,将会获取到相关的信息,如回复“1”的结果如下,

New

至此,公众号建立完毕并成功地与已有的博客关联起来。

欢迎关注晓伟博客的公众号

wxwblog

全国主要机场的航空公司份额计算

近年来航空出行逐渐地普及开来,相应地也产生了许多的问题,比如飞机延误等,其中还牵涉到一些服务态度的问题。这时候有人就说了,航空公司搞垄断、店大欺客等等,其实飞机延误,或者服务态度的问题和垄断并没多大关系。航空业本身就是一个高门槛的烧钱行业,即便国家随便放开行业准入,最后的结果依然会是几家公司独大的局面,因为在这行业里,没有足够的规模是生存不下去的。

说到垄断,对应的就得有判断标准,最直观的标准就是市场份额。比如某某几个公司共同占据了某个行业90%的份额,那这些公司就是垄断的,但至于这些垄断是商业竞争形成的,还是行政管制形成的,则并不在考虑的范围之内。实际上垄断的定义远没有这么简单,否则那些反垄断的官司就不会一个比一个旷日持久了。本文不讨论这个有争议的话题,只单纯研究一下在全国80座主要的机场中,各个航空公司的市场份额问题。

所谓市场份额,必需要有一个标的的指标,比如销售额的市场份额、或者产量的市场份额等等。航空公司通常所涉及的市场份额,一般是指运力的市场份额,即ASK(Available Seat Kilometers),这个指标等于一架飞机上的座位数量乘以其飞行距离所得出的结果。但是由于飞机的实际运营机型并不稳定,因此本文选取了各公司在对应机场的飞行班次作为计算市场份额的指标。

主要的数据来源于PAXIS的运力数据,这里面主要包括飞行班次、座位数以及ASK三个指标,格式如下图所示。第一列为各航空公司的代码,就是我们平时坐飞机时常听到的CA、CZ、MU、HU,图中的3U指的是四川航空公司,CA就是我大擦航了。

PEK_Data

 

 

 

 

由于我们分析的对象是全国主要的80座机场,因此源数据中会包含80个与上表格式类似的Excel表。为了宏程序批处理的实现,这些表统一由机场的三字代码命名,且这80座机场的名单列表也单独存放在一个名为AirportsList的Excel表中,格式如下。

AirportsList

另外,为了结果表的易读性,还需要有航空公司代码和中文名称的对应表,以及不同航空公司对应航空系的列表。所谓航空系指的不同航空公司之间的股权或者经营权隶属关系的集合,比如国航旗下有深航、山航、大连航空,东航则有上海航空、联合航空等等。表的结构如下图所示。

OpAirList

OpAirList2

接下来就码代码实现80座机场中,各航空公司、航空系的班次份额计算,主要的思路如下:

  1. 读入含有80座机场代码和名称的Excel表AirportsList,将数据存入AirList数据集中;
  2. 依次读取AirList数据集中的机场代码,对每座机场作如下运算:
  3. 读入机场对应的班次、座位数、运力数据,(如需输出航空系份额时,在此处对源数据按各航空系进入分类汇总);
  4. 计算各航空公司(航空系)在此机场中的班次份额的座位数份额,计算完成后,继续读取下一个机场的数据进行计算,直接80座机场的数据全部计算完成;
  5. 进行一系列的表合并操作,使得结果更加的直观;

最后还有一个相对市场份额的概念,所谓相对市场份额是指,当某航空公司在一座机场的份额排名为第一时,其相对市场份额便等于他的市场份额与排名第二公司市场份额的比值。当某航空公司在一座机场的份额排名不是第一时,其相对市场份额则等于他的市场份额与排名第一公司市场份额的比值。利用相对市场份额的概念,可以很清晰地将各个公司的优势市场用可视化的方法展示出来。

其中在计算相对市场份额时,遇到了一个小小的插曲:原本打算直接利用Lag函数来计算,但是由于要将数据分为份额排名第一和非排名第一两种情况,因此需要引入if条件语句。然而没有想到的是,lag函数是无法在if这样的条件语句中使用的,具体的原因我没有研究清楚,后来是采取了引入辅助变量的办法解决了这个问题。

具体的计算代码如下,由于初始代码写完后一直未进行修改,所以显得比较凌乱。不过没关系,未来有时间对代码进入修改的时候,修改的过程本身就会变成一个自我提高的过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
*** Macro to calculate airlines group share;
*** dt(parameter) is the source data of each airport;
%macro MergeAirGroup(dt);
    data AirGroup;
    	infile AirShare(OpAirList2.csv) dsd missover firstobs = 2;
    	informat OpAir $2. AirGroup $8.;
    	input OpAir AirGroup;
   run;
   proc sort data = AirGroup; by OpAir; run;
   proc sort data = &dt; by OpAir; run; 
   data &dt;
   	merge &dt AirGroup;
   	by OpAir;
   	if missing(AirGroup) then AirGroup = OpAir;
   run;
   proc sort data = &dt; by AirGroup; run;
   *** Sum the Flights and Seats num for the airlines
   	in each air group;
    	data &dt;
            set &dt(rename = (Flights = Flights2 Seats = Seats2));
            retain Flights Seats;
            by AirGroup;
            if first.AirGroup then do;
            	Flights = 0; Seats = 0;
            end;
            Flights + Flights2; Seats + Seats2;
            if last.AirGroup then output;
            drop Flights2 Seats2;
    	run;
%mend MergeAirGroup;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
*** Macro to get the draft result table;	
%macro AirShare(ns = 3); 
	*** RstDt is the name of the result table;
    	*** ns is the number of obs in each airport result table;
 
    	*** delete and empty the working library;
    	proc datasets lib = work memtype = data kill nolist;
    	run;
    	quit;
    	*** Create a empty dataset to store the result observations;
    	data ShareResult;
	    run;
    	*** Import Aircode matching file;
    	data AirList;
            infile AirShare(AirportsList.csv) dsd missover firstobs=1;
            informat Rank 3. Airport $24. AirCode $3.;
            input Rank Airport AirCode;
    	run;
    	*** Assign the airport names to macro vars Air1 - Air80;
    	data _null_;
            set AirList(keep = AirCode);
            ns = input(_n_, 4.);
            call symput("Air" || compress(_n_), AirCode);
    	run;
 
    	%do i = 1 %to 80 %by 1;
            *** Import the source data of each airport ;
            data &&Air&i;
            	infile AirShare(&&Air&i...csv) dsd missover firstobs = 6;
            	informat MtAir $2. OpAir $2. Flights comma8. Seats comma8.;
            	input MtAir OpAir Flights Seats;
            	if missing(MtAir) then stop;
            run;
    	******************* Merge Group; 
    	%MergeAirGroup(&&Air&i);
    	******************* Merge finished;
            *** Create 2 macro vars to store the sum of the Flights and Seats;
    	    data _null_;
    		        set &&Air&i(keep = Flights Seats) end = eof;
            		SumFlights + Flights; SumSeats + Seats;
            		if eof then do;
                            call symput("SumFlights", SumFlights);
                	    call symput("SumSeats", SumSeats);
            		end;
            run;
	    *** Sort the source table by Flights;
	    proc sort data = &&Air&i; by descending Flights; run;
            *** Calculate the share of Flights and Seats ;
            data &&Air&i..Share;
            	*** Set the order of vars;
            	retain SumFlights SumSeats AirCode OpAir 
	               Flights FlightsShare Seats SeatsShare;       
            	set &&Air&i(keep = OpAir AirGroup Flights Seats obs = &ns);  
            	*parameter ns, control the number of observations for each airport in the result table;
            	if _n_ eq 1 then do;
                    SumFlights = symget("SumFlights");
                    SumSeats = symget("SumSeats");
                    AirCode = "&&Air&i";
            	end;
            	    FlightsShare = Flights / SumFlights;
            	    SeatsShare = Seats / SumSeats;
            	    SubRank = _n_;
            	    format FlightsShare SeatsShare percent8.2;
            	    drop SumFlights SumSeats;
        	run;
        	*** Combine each result table;
        	data ShareResult;
            		set ShareResult &&Air&i..Share;
        	run;
    	%end;
    	*** Sort ShareResult and AirList for Merger;
    	proc sort data = ShareResult; by AirCode; run;
    	proc sort data = AirList; by AirCode; run;
    	*** Merge 2 datasets and delete the missing values;
	    data AirList;
		        merge AirList ShareResult;
        		by AirCode;
        		if missing(Flights) then delete;
    	run;
    	*** Import Operating Air List;
    	data OpAirList;
        		infile AirShare(OpAirList.csv) dsd missover firstobs = 2;
        		informat OpAir $2. Airline $8.;
        		input OpAir Airline;
    	run;
	    	proc sort data = AirList; by OpAir; run;
    	proc sort data = OpAirList; by OpAir; run;
    	data CNRst;
        		retain Rank	Airport	AirCode	OpAir Airline Flights
			            FlightsShare	Seats SeatsShare SubRank;
        		merge OpAirList AirList;
        		by OpAir;
    	run;
    	proc sort data = CNRst; by Rank Aircode descending FlightsShare; run;
    	proc export data = CNRst outfile = "D:\saswork\AirportsShare\CNRst.csv"
                            				dbms = csv replace;
    	run;
%mend AirShare;
*** Get the draft result table without Chinese specification  ;
%AirShare(ns = 10);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
*** Match the chinese name to each airport;
data Airport;
	    set CNRst(keep = Airport);
run;
proc sort data = Airport; by Airport; run;
data Airport;
    	set Airport;
    	by Airport;
    	if first.Airport then output;
run;
 
*** AirportCityRef is an external file with 2 vars Airports and 
	their corresponding chinese names;
data AirportCityRef;
    	infile AirShare(AirportCityRef.csv) dsd missover firstobs = 2;
    	informat Airport $24. City $8.;
	    input Airport City;
run;
 
proc sort data = CNRst; by Airport; run;
proc sort data = AirportCityRef; by Airport; run;
data CityShare;
    	merge AirportCityRef CNRst;
    	by Airport;
    	if missing(city) then delete;
run;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
*** Bostion Index Computing starts;
 
*** File AirportsPaxGrowth has 2 vars, city names and their growth rate
	in recent 3 years ;
data AirportGrowth;
    	infile AirShare(AirportsPaxGrowth.csv) dsd missover firstobs = 2;
    	informat City $8. Growth 8.;
    	input City Growth;
run;
 
*** Merge Airports growth into CityShare by City;
proc sort data = AirportGrowth; by City; run;
proc sort data = CityShare; by City; run;
data CityShare;
    	merge AirportGrowth CityShare;
    	by City;
    	if missing(Rank) then delete;
run;
proc sort data = CityShare; by Rank Aircode descending FlightsShare; run;
********************* Merging finished *****************************;
 
 
*** Create a temp var(bst_temp)  to get the Boston Matrix value;
*** Set all the bst_temp to the biggest flightsvalue of each city;
data CityShare;
    	set CityShare;
    	retain bst_temp;
    if SubRank eq 1 then bst_temp = Flights;
run;
*** Sort cityshare in ascending flights sequence;
proc sort data = CityShare; by Rank Aircode Flights; run;
*** Assign No.2 bst value to No.1's bst value;
data CityShare;
    	set CityShare;
    	bst_temp2 = lag1(Flights);
    	if SubRank eq 1 then bst_temp = bst_temp2;
    	drop bst_temp2;
run;
*** Compute bst index value for each city;
data CityShare;
    	set CityShare;
    	bst = Flights / bst_temp;
    	drop bst_temp;
run;
proc sort data = CityShare; by Rank Aircode descending Flights; run;
proc export data = CityShare outfile = "D:\saswork\AirportsShare\CityShare.csv"
                             			dbms = csv replace;
run;

最后得到整个程序计算出的结果表,

AirportShare

根据这个结果表,就可以画出类似于下图中的相对市场份额的分布表。

BST_Data

利用百度地图计算两城市的距离

最近有个很火的概念叫京津冀经济圈,其中提到北京、河北、天津要做到协同发展。由于工作的关系,经济能够听到协同这个词,同时也知道这是件很难的事情。几家拥有百十来架飞机的公司间要做到协同发展,尚且如此困难,更别提三个省级行政区间的协同了。

京津冀经济圈并不是本文要叙述的重点,本文只对“圈”字感兴趣。实际上这个概念并不陌生,什么珠三角经济圈、长三角经济圈、环渤海经济圈等等,都是在围绕一个圈在做文章。那么这个圈应当如何界定呢?恰好在前不久的工作中,我接触到了一个机场辐射圈的问题。大致意思是以机场为圆心,在一定公里数范围内的地区,就可以认为是这个机场可能影响到的地区。其实简化来讲,就是要解决某个机场,比如首都机场,在它的一定公里半径范围内,都覆盖到了哪些城市和地区。定量地说,就是要批量计算不同城市对间距离的问题。

确定了问题之后,我开始查找相关的资料,之后便遇到了一个有意思的问题。我们经常听到这样的说法,比如山西省长治市与北京市的直线距离为550公里。咋一听这个说法好像没什么问题,但是仔细一想,所谓距离应该指的是两点间的距离,但是北京和长治显然不是两个点,长治市到房山的距离和到怀柔的距离肯定不一样,那么这个550公里到底是哪两个点间的距离呢?这个问题提醒我要计算城市间的距离,必然不可能做到非常精确,我所要做的应该是找到每个城市的“地标”点,然后计算这些点间的距离。但是城市的地标点怎么确定呢?

北京的地标自然是好找的,但是某些不太知名的城市,就不是那么容易了。思来想去,突然想到一句话,普天之下,莫非王土;率土之滨,莫非王臣。小的城市肯定没有天安门, 也没有奥林匹克公园,但是每个城市肯定有个叫某某市人民政府的地方,所以地标点就确定为每个城市人民政府的所在地。

之后就是另一个关键性的问题,怎么样完成从坐标到距离的换算呢?由于坐标需要从百度地图之类的网站上截取下来的,因此就需要寻找由地图坐标到实际距离的换算公式,在网上搜寻一番后发现还真有相关的计算公式,于是操起SAS,开始流程化作业,流水账如下:

1,登入百度地图API开发页面中的坐标拾取工具,逐个搜索每一个地标点的坐标,然后复制到源数据文件中。此处方法显得粗暴且笨拙,我猜想大概有抓取数据文件的批处理程序,但是没有深入去研究。

CoordinatePic.png

 

 

 

 

一番机械的重复劳动后,得到了全国289座城市的市政府坐标汇总数据表,同样的方法可以得到全国80座主要机场的坐标数据。

CityGovernment.png

 

 

 

 

注意此时的坐标是作为一个变量保存下来的,后面的程序处理中,要从这个变中提取出X轴坐标和Y轴坐标。

2,码程序。大致的思路如下:

  1. 读入相关源数据(废话),包括机场、城市的坐标数据,以及各城市的GDP、人口数量以及人均收入数据;
  2. 将GDP、人口数量、人均收入三个数据集,以城市名称为By变量,合并入同一表中(CityInfo)。这个表中储存了289个城市的名称,以及对应的经济指标数据;
  3. 依靠宏变量(Nair),对80座机场进行遍历,将每一座机场分别与289个城市数据合并,并放入临时数据集tempobs中,最终将所有机场数据放入数据集AirportRegion中,AirportRegion大约有80*289条观测;
  4. 对AirportRegion中的每条观测应用坐标到距离的换算公式,筛选掉半径范围(Radius)之外的数据;

具体的代码如下:第二段程序的第40-48行就是计算公式的表达式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
********************************************************************;
*                  Macro Name:     CityRadius
*                  Author:         Wang Xiaowei
*                  Function:       Get Regional GDP Population for
*                                a given radius of an airport
*                CreateDate:       07Feb2014
********************************************************************;
%include "d:/saswork/ShareFunctions.sas";
 
/* Work space path */
filename Radius "d:/saswork/Radius/";
 
/* Source data preparation */
options pageno = 1 linesize = 100;
/* GDP data for each city */
data CityGDP;
    infile  Radius(CityGDP.csv) dsd missover firstobs = 2;
    informat Province $6. Unit $10. City $10.
            YC2009 - YC2012 8. Growth percent8.;
    input Province Unit City YC2009 - YC2012 Growth;
run;
/* Population data for each city */
data CityPopulation;
    infile Radius(CityPopulation.csv) dsd missover firstobs = 2;
    informat Province $6. City $10. Unit $4.
            CityX009 - CityX012 8. Growth percent8.;
    input Province City Unit CityX009 - CityX012 Growth;
run;
/* Disposable income data for each city */
data CityIncome;
    infile Radius(CityIncome.csv) dsd missover firstobs = 2;
    informat Province $6. Unit $6. City $10.
            YC2009 - YC2012 8. Growth percent8.;
    input Province Unit City YC2009 - YC2012 Growth;
run;
/* Coordinate data for each city */
data AirportCoordinate;
    infile Radius(AirportCoordinate.csv) dsd missover firstobs = 2;
    informat Airport $20. City $10. AirportLC $25.;
    input Airport City AirportLC;
    /*    Extract the X and Y num from var AirportLC */
    AirportX = put(scan(compress(AirportLC), 1, ","), 10.);
    AirportY = put(scan(compress(AirportLC), 2, ","), 10.);
    drop AirportLC;
run;
*proc contents data = AirportCoordinate;
*run;
data CityCoordinate;
    infile Radius(CityCoordinate.csv) dsd missover firstobs = 2;
    informat Province $6. City $10. CityLC $25.;
    input Province City CityLC;
    /* Extract the X and Y num from var AirportLC */
    CityX = put(scan(compress(CityLC), 1, ","), 10.);
    CityY = put(scan(compress(CityLC), 2, ","), 10.);
    drop CityLC;
run;
 
/* Sort 3 datasets to merge them together */
%sort(CityGDP, City);
%sort(CityPopulation, City);
%sort(CityIncome, City);
%sort(CityCoordinate, City);
data CityInfo;
    merge   CityGDP(drop = YC2009 - YC2011 Growth rename = (YC2012 = GDP2012))
            CityPopulation(drop = CityX009 - CityX011 Growth
            rename = (CityX012 = Population2012))
            CityIncome(drop = YC2009 - YC2011 Growth
            rename = (YC2012 = Income2012))
            CityCoordinate;
    by City;
    drop Unit;
run;
/* Source Data has been prepared successfully*/
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
/* Radius(parameter) is the distance between the airport and city */
/* The unit of radius is meter, not km */
%macro DistCal(radius = 200000);
    /* Create an empty dataset to store the results */
    data AirportRegion;
    run;
    /* Nairport(macro var) has the number of airports
        in the dataset AirportCoordinate */
    /* Ncity(macro var) similar to nair */
    data _null_;
        if 0 then do;
            set AirportCoordinate nobs = nair;
            set CityInfo nobs = ncity;
        end;
        call symput("Nairport", nair);
        call symput("Ncity", ncity);
    run;
    %do sair = 1 %to &Nairport;
        /* Extract the &sair(th) airport */
        data AirportTemp;
            set AirportCoordinate;
            if _n_ eq &sair;
        run;
        /* Repeate the &sair(th) airport data for &Ncity(289) times */
        data AirportTemp2;
            set %do i = 1 %to &Ncity; AirportTemp %end;;
        run;
        /* Combine the repeated airport data to the city data */
        data tempobs;
            merge AirportTemp2 CityInfo(rename = (City = Cit2));
        run;
        /* Store the temp data to the result dataset(AirportRegion) */
        data AirportRegion;
            set AirportRegion tempobs;
        run;
        %end; /* sair end */
     /* Calculation procedure  */
    data AirportRegion;
        set AirportRegion;
            R = 6370996.81; /* R is the radius of earth */
            pk = 180 / constant("pi");
            a1 = AirportY / pk; a2 = AirportX / pk;
            b1 = CityY / pk; b2 = CityX / pk;
            /*        Calculation formula */
            t1 = cos(a1) * cos(a2) * cos(b1) * cos(b2);
            t2 = cos(a1) * sin(a2) * cos(b1) * sin(b2);
            t3 = sin(a1) * sin(b1);
            dist = R * arcos(t1 + t2 + t3) * 1.24;
         /* 1.24 is a ratio of straight dist and real dist*/
        drop R pk a1 a2 b1 b2 t1 t2 t3
            AirportX AirportY CityX CityY;
        /* Drop the cities outside the radius area */
        if _n_ ge 1 and dist le &radius;
        dist = round(dist / 1000, 1);
    run;
%mend DistCal;
%DistCal(radius = 200000);
proc export data = AirportRegion outfile = "d:/saswork/Radius/AirportRegion.csv"
                    dbms = csv replace;
run;
/* End */

最后得到数据结果,以200KM范围为例

AirportsArea.png

 

 

 

 

最后要说明的是,这个计算距离,与实际的驾车距离间仍然有所偏差,二者的比例大约为1.24左右,当然这个比例可能会因城市的不同而异,毕竟每个城市的实际交通状况是有差异的。