大数据分析案例
数据分析案例
——发掘数据背后的真实故事，发现规律，发现那些意想不到的规律.......
（更新，）
实际案例分类：
1.1 金融案例1：已知某种保险产品在一个保单年度内损失情况（数据如下表所示），其中给出了不同损失次数下的保单数，现对损失次数的整体规律进行分析。
num=c(rep(0:5,c(,41,10,4)))
lambda=mean(num)
ppois=dpois(k,lambda)
poisnum=ppois*length(num)
plot(k,poisnum,ylim=c(0,1600))
samplenum=as.vector(table(num))
points(k,samplenum,type=&p&,col=2)
legend(4,1000,legend=c(&num&,&poisson&),col=1:2,pch=&o&)
本案例服从参数Lambda=0.4780571的泊松分布，可以此开展相关预测。
1.2 金融案例2：某银行规定VIP客户的月均账户余额要达到100万元，并以此作为比较各分行业绩的一项指标，从三个分行中，分别随机抽取7个VIP客户的账户余额，首先对三个分行的账户余额分别进行正态检验。
& a1=c(103,101,98,110,105,100,106)
& a2=c(113,107,108,116,114,110,115)
& a3=c(82,92,84,86,84,90,88)
& shapiro.test(a1)
Shapiro-Wilk normality test
W = 0.9778, p-value = 0.948& shapiro.test(a2)
Shapiro-Wilk normality test
W = 0.9189, p-value = 0.4607& shapiro.test(a3)
Shapiro-Wilk normality test
W = 0.9547, p-value = 0.7724
利用Shapiro-Wilk正态检验方法，对三个分行的数据进行分析，结果表明三个P值均大于显著性水平alpha=0.05，说明三组数据都服从正态分布规律，可以作为业绩考核参考指标。
1.3 金融案例3：以国内10家股份制商业银行为研究对象，选取11种指标来研究它们的经营绩效，数据来自《中国金融统计年鉴》，具体数值见下表，试找出反映银行经营绩效的指标体系。
负债资产率
不良贷款率
贷款呆账准备率
资产利润率
资金周转率
投资收益率
营业收入费用率
营业收入增长率
营业支出增长率
利润增长率
& factor.analysis=function(x,m)
+ p=nrow(x)
+ x.diag=diag(x)
+ sum.rank=sum(x.diag)
+ rowname=paste(&X&,1:p,sep=&&)
+ colname=paste(&Factor&,1:m,sep=&&)
+ A=matrix(0,nrow=p,ncol=m,dimnames=list(rowname,colname))
+ eig=eigen(x)
+ for(i in 1:m)
+ A[,i]=sqrt(eig$values[i])*eig$vectors[,i]
+ var.A=diag(A%*%t(A))
+ rowname1=c(&SS loadings&,&ProportionVar&,&Cumulative Var&)
+ result=matrix(0,nrow=3,ncol=m,dimnames=list(rowname1,colname))
+ for(i in 1:m)
+ result[1,i]=sum(A[,i]^2)
+ result[2,i]=result[1,i]/sum.rank
+ result[3,i]=sum(result[1,1:i])/sum.rank
+ method=c(&Principal Component Method&)
+ list(method=method,loadings=A,var=cbind(common=var.A,specific=x.diag-var.A),result=result)
& bank=read.table(&d:/R311/data_23714/bank.txt&,header=T)
& bank=bank[,-1]
& bank_cor=cor(bank)
& options(digits=3)
& factor.analysis(bank_cor,5)
[1] &Principal Component Method&
Factor1 Factor2 Factor3
Factor4 Factor5
-0.209 -0.8
0.911 -0.6 -0.28
0.200 -0.2
0.856 -0.7 -0.10
common specific
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings
0.9424
ProportionVar
0.0857
Cumulative Var
0.9454
使用自己编写的函数factor.analysis()，通过主成分分析方法，将上述11个指标综合成为5个因子，来概括反映银行经营绩效的指标体系。从Cumulative Var可知，这5个因子的方差贡献率累计达到94.54%，因此可以认为这5个因子已经包含了总体的大部分信息。5个主成分（因子）和原始变量的线性关系如下：
F1=0.394*X1-0.209*X2+0.804*X3+0.877*X4+0.594*X5+0.309*X6+0.911*X7+0.200*X8+0.856*X9-0.764*X10+0.644*X11
依此类推，得到F2~F5的表达式。
F1与X3（不良贷款率）、X4（贷款呆账准备率）、X7（投资收益率）、X9（营业收入增长率）呈现高度相关（系数均大于0.8），因此F1反映了银行成长能力；
F2与X6（资金周转率）、X8（营业收入费用率）相关性较高，因此F2反映了银行的经营效率；
F3与X2（负债资产率）、X11（利润增长率）相关性较高，因此F3反映了银行的盈利水平；
F4与X1（流动比率）相关性较高，因此F4反映了银行的流动性和短期偿债能力；
F5与X5（资产利润率）相关性较高，因此F5反映了银行的资金使用情况。
1.4 金融案例4：已知我国2002年1月至2013年8月的货币供应量M1的数据，见下表，试预测其未来发展趋势。
& data=read.table(&d:/R311/data_23714/M1.txt&,header=T)
& M1=ts(data$M1,frequency=12,start=c(2002,1))
& plot.ts(M1)
& library(TTR)
& M1.SMA3=SMA(M1,n=3)& plot(M1.SMA3)
---------------------------------------------------------------------------------
& data=read.table(&d:/R311/data_23714/M1.txt&,header=T)
& M1_2=ts(data$M1,frequency=12,start=c(2002,1))
& M1_2.pre=HoltWinters(M1,gamma=FALSE)
& M1_2.pre
Holt-Winters exponential smoothing with trend and without seasonal component.
HoltWinters(x = M1, gamma = FALSE)
Smoothing parameters:
alpha: 0.6630848
beta : 0.125953
gamma: FALSE
Coefficients:
& plot(M1_2.pre)
& library(forecast)
& M1_2.pre2=forecast.HoltWinters(M1_2.pre,h=20)
& plot.forecast(M1_2.pre2,main=&用HoltWinters模型对将来5年（20个季度）的M1值进行预测&)
& acf(M1_2.pre2$residuals,lag.max=10)
& Box.test(M1_2.pre2$residuals,lag=10,type=&Ljung-Box&)
Box-Ljung test
M1_2.pre2$residuals
X-squared = 22.1802, df = 10, p-value = 0.01421
-------------------------------------------------------------------------------------
& data=read.table(&d:/R311/data_23714/M1.txt&,header=T)
& M1_3=ts(data$M1,frequency=12,start=c(2002,1))
& M1_3.diff=diff(M1_3)
#进行一阶差分
& plot(M1_3.diff,main=&对M1进行一阶差分后的数据曲线&)
-------------------------------------------------------------------------------------
& logM1.diff=diff(log(M1_3))
& op=par(mfrow=c(2,1))
& plot(logM1.diff,main=&先对M1取对数（消除一定的异方差性）再差分效果图&)
& acf(logM1.diff,lag.max=10)
------------------------------------------------------------------------------------
& auto.arima(log(M1_3),ic=&bic&) #ARIMA可以含有差分，故只需事前进行取对数即可
Series: log(M1_3)
ARIMA(0,1,0)(1,0,1)[12]
Coefficients:
sigma^2 estimated as 0.0002146:
log likelihood=380.56
AIC=-755.12
AICc=-754.94
BIC=-746.32
& logM1.arima=arima(log(M1_3),order=c(1,1,1))
& logM1.arima
Series: log(M1_3)
ARIMA(1,1,1)
Coefficients:
sigma^2 estimated as 0.0004199:
log likelihood=341.74
AIC=-677.49
AICc=-677.31
BIC=-668.68
-----------------------------------------------------------------------------
& logM1.pre3=forecast.Arima(logM1.arima,h=12)
& logM1.pre3
& plot.forecast(logM1.pre3)
----------------------------------------------------------------------------
& acf(logM1.pre3$residuals,lag.max=5)
& Box.test(logM1.pre3$residuals,lag=5,type=&Ljung-Box&)
Box-Ljung test
logM1.pre3$residuals
X-squared = 10.6986, df = 5, p-value = 0.05769
----------------------------------------------------------------------------
& M1.pre3=exp(logM1.pre3$mean)
a. 从“Time-M1”示意图可以看出，我国货币供应量M1具有明显的上升趋势，但不存在明显的季节性变动，因此初步判断这个时间序列可以被描述为一个加法模型，使用简单移动平均平滑来估计趋势部分，尝试了跨度为3的平滑，如图&Time-M1.SMA3&，跨度越大越可以消除不规则波动，简单移动平均法展现的趋势部分就越平滑，发展趋势看起来更加清晰。
b. 目测估计M1可能属于平稳性时间序列，选择简单指数平滑法HoltWinters模型对M1进行短期预测，从图形效果来看，原始数据和拟合（预测）值差别不大，目测感觉效果不错。对未来五年（20个季度）的M1值进行预测，预测图中的浅灰色阴影为置信水平90%的预测区间，深灰色阴影为置信水平80%的预测区间，中间粗线条为未来5年的M1预测水平（估计值）。这些拟合（预测）效果是否真的不错呢？下面用残差的自相关图来检验。从残差的自相关图中可以看出，第6阶和第8阶残差自相关系数稍微超过了置信界限，说明残差序列具有一定的相关性。接下来再用Ljung-Box检验来做最终的判断，而Ljung-Box检验结果中p-value = 0.01421，不大于0.05显著性水平，可以拒绝原假设，说明残差序列存在自相关性，用HoltWinters模型拟合（预测）的时间序列信息并不完全。
c. 选择适用于非平稳时间序列的ARIMA（差分自回归移动平均）模型对M1进行拟合和预测。首先对M1进行一阶差分，从差分后的曲线图中可以看出，一阶差分去掉了趋势，变得更平稳了，但是曲线尾端（最后几年）波动性较大，意味着存在“异方差性”。
于是改为先对M1取对数（消除一定的异方差性），然后再进行一阶差分，此时数据变得更加平稳，尾端波动也不再剧烈，对其进行自相关图检验，序列在10阶内自相关系数基本上都落在置信区间内，可以确认其平稳性（注意：这是对数后再差分结果的平稳性）。
e. 平稳化完成之后，根据平稳序列的“偏相关系数和自相关系数”特点来选择合适的模型:
--------(1)若平稳序列的偏相关系数是截尾的，而自相关系数是拖尾的，则适合AR模型；
--------(2)若平稳序列的偏相关系数是拖尾的，而自相关系数是截尾的，则适合MA模型；
--------(3)若
平稳序列的偏相关系数和自相关系数都是拖尾的，则适合ARMA模型；
--------(4)可以用pacf(diff(log(M1)))来观察偏相关系数、acf(diff(log(M1)))观察自相关系数，也可直接用auto.arima(log(M1))自动探寻合适的ARIMA模型，只是耗时较长。注意：pacf()和acf()处理的是diff(log(M1))，而auto.arima()处理的是log(M1)，因为arima自身可包含差分操作。
本案例中，auto.arima(log(M1))自动选出了ARIMA(0,1,0)(1,0,1)模型p=1,d=1,q=1。
f. 利用e.中得到的模型参数，对log(M1)进行预测，
可以看到其具体预测值及其80%、95%的置信区间，也可以看到预测结果曲线，曲线的浅灰色阴影代表90%置信水平的预测区间，深灰色阴影代表80%置信水平的预测区间，中间粗线代表预测水平（具体数值）。
g. 最后，验证ARIMA(1,1,1)预测的有效性，通过acf()的残差自相关图检验和Ljung-Box检验，可知残差的自相关系数均落入置信界限内，初步说明残差不相关；Ljung-Box检验的p-value = 0.05769大于0.05显著性水平，再次确认残差无相关性，属于白噪声，ARIMA(1,1,1)适用于对log(M1)的拟合和预测。
h. 注意，f.、g.步骤中，只是对log(M1)的预测，若要得到真正的对M1的预测，只需将log(M1)的预测值或置信区间，通过exp()函数取对数即可，见上文最后一个exp()函数示意图。
2.1 喷泉案例：某地质学家记录了美国黄石公园内Old Faithful喷泉在一年内的299组喷发数据，现对喷泉整体喷发规律进行分析。
喷发相隔时间
喷发持续时间
library(MASS)
attach(geyser)
ll=function(para)
f1=dnorm(waiting,para[2],para[3])
f2=dnorm(waiting,para[4],para[5])
f=para[1]*f1+(1-para[1])*f2
tmp=sum(log(f))
return(-tmp)
geyser.est=nlminb(c(0.5,50,10,80,10),ll,lower=c(0.0001,-Inf,0.0001,-Inf,0.0001),upper=c(0.9999,Inf,Inf,Inf,Inf))
options(digits=3)
geyser.est$par
p=geyser.est$par[1]
mu1=geyser.est$par[2];sigma1=geyser.est$par[3]
mu2=geyser.est$par[4];sigma2=geyser.est$par[5]
x=seq(40,120)
f=p*dnorm(x,mu1,sigma1)+(1-p)*dnorm(x,mu2,sigma2)
hist(waiting,freq=F)
lines(x,f)
本案例服从具有双峰的两正态混合分布，f(x)=p.N(x;mu1,sigma1)+(1-p)*N(x;mu2,sigma2)，其中参数：p=0.308,mu1=54.203,sigma1=4.952,mu2=80.360,sigma2=7.508
3.1 零售案例1：百货公司Bamberger's是一家为社区提供大众商品的零售商店，长期以来在顾客中有着服务极好的口碑，为了努力维持商店的良好声誉，公司实施了一项新计划——将营业时间延长至深夜。以计划实施前后27周销售额数据为依据，分析评估本计划的实施效果。
计划实施前
计划实施后
sales=read.table(&d:/R311/xxf_R/xxf_data/sales.txt&,header=T)
attach(sales)
par(mfrow=c(1,2))
hist(prior)
hist(post)
twosample.ci=function(x,y,alpha,sigma1,sigma2)
n1=length(x)
n2=length(y)
xbar=mean(x)-mean(y)
z=qnorm(1-alpha/2)*sqrt(sigma1^2/n1+sigma2^2/n2)
c(xbar-z,xbar+z)
twosample.ci(post,prior,alpha=0.05,sd(prior),sd(post))
var.test(prior,post)
新计划实施后，Bamberger's公司的周营业额明显增加，增加幅度的范围是[18.80,30.28]；同时周营业额的波动性变大[1.17, 5.64]倍。
3.2 零售案例2：某种饮料有五种不同的包装方式，分别在五个不同地区销售，现从每个地区随机抽取一个规模相同的超市，得到该饮料不同包装的销售数据，见下表，试分析不同包装和不同地区对销售额的影响。
& x=c(20,12,20,10,14,22,10,20,12,6,24,14,18,18,10,16,4,8,6,18,26,22,16,20,10)
& sales=data.frame(x,A=gl(5,5),B=gl(5,1,25))
& bartlett.test(x~A,data=sales)
Bartlett test of homogeneity of variances
Bartlett's K-squared = 0.6653, df = 4, p-value = 0.9555& bartlett.test(x~B,data=sales)
Bartlett test of homogeneity of variances
Bartlett's K-squared = 1.2046, df = 4, p-value = 0.8773
& sales.aov=aov(x~A+B,data=sales)
& summary(sales.aov)
Df Sum Sq Mean Sq F value Pr(&F)
2.303 0.1032
3.874 0.0219 *
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Bartlett检验结果中，因素A和因素B的P值都远大于0.05显著性水平，说明因素A、B的样本数据满足方差齐性，可以继续进行双因素方差分析；
双因素方差分析结果中，因素A的P=0.，因素B的P=0.，说明包装对销售的影响不明显，地区对销售有显著影响。
4.1 市政实例1：某市为了了解居民住房情况，抽查了2000户家庭，其中人均不足5平米的困难户有214个，依据以上数据分析该市困难户比率的大致范围。
n=2000,x=214
binom.test(214,2000)
Exact binomial test
214 and 2000
number of successes = 214, number of trials = 2000, p-value&
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
sample estimates:
probability of success
城市人口数较大，本案例抽样比率较小，故采用二项式检验，结果表明困难户比率的95%置信水平下的区间范围为[9.38%,12.14%]。
4.2 市政实例2：某市2012年各月新建住宅价格指数已知（见下表），要检验此指数是否服从均值为102.4的正态分布规律，并检验方差是否不超过0.25。
&bj=c(102.5,102.4,102.0,101.8,101.8,102.1,102.3,102.5,102.6,102.8,103.4,104.2)
t.test(x=bj,mu=102.4)
One Sample t-test
t = 0.6701, df = 11, p-value = 0.5166
alternative hypothesis: true mean is not equal to 102.4
95 percent confidence interval:
sample estimates:
&chisq.var.test=function(x,var,mu=Inf,alternative=&two.sided&)
n=length(x)
if(mu&Inf){df=n;v=sum((x-mu)^2)/n}
chi2=df*v/var
options(digits=4)
result=list()
result$df=result$var=v;result$chi2=chi2;
result$P=2*min(pchisq(chi2,df),pchisq(chi2,df,lower.tail=F))
if(alternative==&greater&) result$P=pchisq(chi2,df,lower.tail=F)
else if(alternative==&less&) result$P=pchisq(chi2,df)
chisq.var.test(bj,0.25,alternative=&less&)
t统计量检验结果P=0.5166&alpha=0.05，因此认为2012年该市各月新建住宅价格指数服从均值为102.4的正态分布；
卡方检验的结果P=0.9656远大于alpha=0.05，说明新建住宅价格指数的方差大于0.25，即价格指数变动很大。
4.3 市政实例3：根据指标对全国的主要城市进行分类，从而帮助政府有针对性地对不同城市制定政策。下表是2012年26座城市的5项经济指标，分别是：平均房价、人均GDP、平均工资、人均可支配收入和CPI指数，试对这些城市进行分类，以便制定不同的政策。
人均可支配收入
& dat=read.table(&d:/R311/data_23714/real_estate.txt&,header=T)
& dat=dat[,-1]
& d=dist(dat)
& hc=hclust(d,method=&ward.D&)
& plot(hc,hang=-1)
利用“离差平方和法”聚类分析，根据房价和经济指标的相似性，将26座城市分成了四类，如示意图所示。第一组城市的特点是人均工资水平都不算太高，而房价也处于较低水平；类似的，对其它各组也能找出共同点，以此为基础，国家宏观调控部门可以给出不同城市应当采取的房价政策。
5.1 医学案例：为研究某药物对高血压病人的治疗效果，用20名患者分成两组开展试验，试问该药物是否有降压作用。
待试药物组
& x=c(117,127,141,107,110,114,115,138,127,122)
& y=c(113,108,120,107,104,98,102,132,120,114)
& t.test(x,y,paired=TRUE,alternative=&greater&)
Paired t-test
t = 4.586, df = 9, p-value = 0.0006586
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
sample estimates:
mean of the differences
t检验结果显示P=0.0006586，远小于alpha=0.05，说明待试药物组的血压均值明显降低，即该药物有降压作用。
6.1 工业实例：对某台设备进行寿命检验，记录10次无故障工作时间（见下表），检验其是否服从参数为1/1500的指数分布规律。
& x=c(420,500,920,50,00,2350)
& ks.test(x,&pexp&,1/1500)
One-sample Kolmogorov-Smirnov test
D = 0.3015, p-value = 0.2654
alternative hypothesis: two-sided
KS检验结果P=0.2654，远大于alpha=0.05，说明该设备寿命服从参数为1/1500的指数分布。
7.1 教育实例1：某制造商雇佣了来自三所本地大学的毕业生作为管理人员，公司人事部门收集信息并考核了年度工作成绩，从三所大学来的职员中随机抽取了独立样本，样本量分别为7、6、7，数据如下表所示，制造商想知道来自这三所大学的毕业生在管理岗位上的表现是否有所不同。
& data=data.frame(x=c(25,70,60,85,95,90,80,60,20,30,15,40,35,50,70,60,80,90,70,75),g=factor(rep(1:3,c(7,6,7))))
& kruskal.test(x~g,data=data)
Kruskal-Wallis rank sum test
Kruskal-Wallis chi-squared = 8.9839, df = 2, p-value = 0.0112
Kruskal-Wallis秩和检验的结果P=0.0112&alpha=0.05，说明来自这三个大学的毕业生在管理岗位上的表现有比较显著的差异。
7.2 教育实例2：一项研究汉字读写能力与数学的关系的调查，获取了232个美国亚裔学生的数学成绩和汉字读写能力的数据，见下表，其中汉字读写能力指标有3个水平：“纯汉字”表示可以流利自如地进行纯汉字读写，“半汉字”表示读写中只有部分汉字，而“纯英文”表示只能读写英文而不懂汉字；数学成绩用A、B、C、D四个水平来表示，试分析汉语读写能力与数学成绩有何关联性。
汉字读写能力
& ch=data.frame(A=c(47,22,10),B=c(31,32,11),C=c(2,21,25),D=c(1,10,20))
& rownames(ch)=c(&Pure-Chinese&,&Semi-Chinese&,&Pure-English&)
& library(MASS)
& ch.ca=corresp(ch,nf=2)
& options(digits=4)
First canonical correlation(s): 0.9
Row scores:
Pure-Chinese
Semi-Chinese -0.9
Pure-English -1.0
Column scores:
& biplot(ch.ca,main=&汉语读写与数学成绩之间关联性分析示意图&)
关联分析（对应分析）是一种可视化的多元统计方法，主要通过图形分析来得出结论，从“汉语读写与数学成绩之间关联性分析示意图”中两种散点的横坐标之间的距离（纵坐标的距离对于分析意义不大）可以看出“纯汉字”和数学成绩A最接近，说明数学好的人大多可以自如地进行纯汉字读写，或者说可以进纯汉字读写的人大多数学成绩较好；“纯英文”和数学成绩D非常接近，说明数学差的人大多不会汉字只会英文，或者说不会汉字只会英文的人大多数学成绩较差；而“半汉字”介于数学成绩B和C之间，说明会部分汉字的学生，其数学成绩一般。
8.1 交通实例：某城市道路交通管理部门为了研究不同路段、不同时段的拥堵情况，分别在三个路段和高峰期、非高峰期进行试验，每一个水平组合下测量5次，共获得30个行车时间的数据，如下表所示，试分析路段和时段对行车时间的影响。
25 24 27 25 25
19 20 23 22 21
29 28 31 28 30
20 17 22 21 17
18 17 13 16 12
22 18 24 21 22
& data0=c(25,24,27,25,25,19,20,23,22,21,29,28,31,28,30,20,17,22,21,17,18,17,13,16,12,22,18,24,21,22)
& traffic=data.frame(data0,A=gl(2,15,30),B=gl(3,5,30,labels=c(&I&,&II&,&III&)))
& bartlett.test(data0~A,data=traffic)
Bartlett test of homogeneity of variances
data0 by A
Bartlett's K-squared = 0.0533, df = 1, p-value = 0.8174& bartlett.test(data0~B,data=traffic)
Bartlett test of homogeneity of variances
data0 by B
Bartlett's K-squared = 0.5776, df = 2, p-value = 0.7492& op=par(mfrow=c(1,2))
& plot(data0~A+B,data=traffic)
& attach(traffic)
& interaction.plot(A,B,data0,legend=F)
& interaction.plot(B,A,data0,legend=F)
& result.aov=aov(data0~A*B,data=traffic)
& summary(result.aov)
Df Sum Sq Mean Sq F value
1 313.63
313.63
84.766 2.41e-09 ***
2 261.60
130.80
35.351 7.02e-08 ***
3.33
0.42
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
（1）Bartlett检验结果中时段因素A和路段因素B的相应P值均远大于显著性水平0.05，说明两个因素的数据满足方差齐性要求，可以进一步做方差分析；
（2）箱线图表明，单独观察时段因素和路段因素对行车时间的影响，发现因素的不同水平对行车时间都有明显影响；
（3）从interaction.plot()所绘制的时段和路段交互效应图A、B可初步判断：时段和路段两种因素之间没有交互效应；
（4）利用方差分析，对（2）、（3）中图形判断结果进行确认，因素A和因素B的P值均远远小于显著性水平0.05，而因素交互A:B的P值大于0.05，说明时段和路段对行车时间有显著影响，但两个因素之间没有明显的交互效应。
8.2 交通案例2：对于某车险保单，考虑车型、被保险人性别等因素对索赔次数的影响。
风险暴露数
& dat=data.frame(y=c(42,37,10,101,73,14),n=c(500,,500,300),type=rep(c('小','中','大'),2),gender=rep(c('男','女'),each=3))
& dat$logn=log(dat$n)
& dat.glm=glm(y~type+gender,offset=logn,data=dat,family=poisson(link=log))
& summary(dat.glm)
glm(formula = y ~ type + gender, family = poisson(link = log),
data = dat, offset = logn)
Deviance Residuals:
Coefficients:
Estimate Std. Error z value Pr(&|z|)
(Intercept)
& 2e-16 ***
1.1e-08 ***
& 2e-16 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 163.490
degrees of freedom
Residual deviance:
degrees of freedom
AIC: 61.677
Number of Fisher Scoring iterations: 5
从分析结果P值可知，小车、女士对应的P值远远小于0.05，说明显著，而中型车虽然小于0.05但不是远远小于，说明比较显著，即小车、女士的索赔次数非常显著，而中型车的索赔次数比较显著。
9.1 农林案例1：为研究三种肥料对苹果树产量的影响有无明显差异，将每种肥料各施用于8棵果树，测量每棵树的初始产量和施肥后增量，如下表所示，试进行分析判断。
果树和肥料
15,13,11,12,12,16,14,17
17,16,18,18,21,22,19,18
22,24,20,23,25,27,30,32
A组果树+a种肥料
85,83,65,76,80,91,84,90
B组果树+b种肥料
97,90,100,95,103,106,99,94
C组果树+c种肥料
89,91,83,95,100,102,105,110
Weight_Initial=c(15,13,11,12,12,16,14,17,17,16,18,18,21,22,19,18,22,24,20,23,25,27,30,32)
Weight_Increment=c(85,83,65,76,80,91,84,90,97,90,100,95,103,106,99,94,89,91,83,95,100,102,105,110)
feed=gl(3,8,24)
data_feed=data.frame(Weight_Initial,Weight_Increment,feed)
library(HH)
m=ancova(Weight_Increment~Weight_Initial+feed,data=data_feed)
Analysis of Variance Table
Response: Weight_Increment
Sum Sq Mean Sq F value
Pr(&F)
Weight_Initial
1 1621.12 1621.12 142.445 1.496e-10 ***
31.071 7.322e-07 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
协方差分析结果的P值远远小于0.05，说明试验结果有显著差别，即三种肥料对苹果产量有很大的影响；从图中纵向位置也可以看出，在三种不同肥料作用下，苹果产量增量有所不同，B、C肥料的苹果增量明显高于A肥料；从图中横向位置也可以发现，C苹果树的初始产量（施肥之前产量）本事就高于A果树。
农林案例2：某地区农业生态经济系统中，各子区域单元相关经济指标数据如下表所示，请将数据进行综合，以较更典型的指标信息来描述该地区农业生态经济的发展状况。
人口密度（人/km^2）
人均耕地面积（ha）
森林覆盖率（%）
农民人均纯收入（元/人）
人均粮食产量（kg/人）
经济作物占农作物播种面积比例（%）
耕地占土地面积比率（%）
果园与林地面积之比（%）
灌溉田占耕地面积之比（%）
& agri=read.table(&d:/R311/data_23714/agriculture.txt&,header=TRUE)
& agri=agri[,-1]
& agri.pr=princomp(agri,cor=TRUE)
# cor=TRUE:利用相关阵计算
& options(digits=4)
& summary(agri.pr,loadings=TRUE)
& screeplot(agri.pr,type=&line&,main=&主成分示意图&)
& biplot(agri.pr,choices=1:2,main=&主成分散点图&)
& biplot(agri.pr,choices=2:3,main=&主成分散点图&)
Importance of components:
Comp.1 Comp.2 Comp.3
Comp.5 Comp.6
Standard deviation
2.5 1.34 0.9 0.33821
Proportion of Variance 0.2 0.38 0.4 0.01271
Cumulative Proportion
0.1 0.34 0.8 0.99147
Standard deviation
Proportion of Variance 0..003497
Cumulative Proportion
X1 X2 X3 X4 X5 X6 X7 X8 X9
0.342 -0.368
＿ -0.375 -0.355
＿ 0.614
＿ 0.155 -0.761 -0.110 ＿ ＿ ＿
＿ ＿ ＿ ＿ 0.206
0.467 -0.203 -0.692
＿ 0.601
＿ -0.598
＿ ＿ 0.139
＿ ＿ 0.396 -0.508
0.580 ＿ ＿
0.638 ＿ ＿ ＿
＿ -0.246 -0.148
＿ ＿ -0.241 -0.777 -0.235
＿ ＿ 0.950
＿ ＿ ＿ ＿ -0.231 ＿
＿ ＿ -0.224
＿ -0.136 -0.246
0.532 -0.613
从方差贡献率Proportion of Variance数据可以看出，第一主成分Comp.1的贡献率为51.79%，第二主成分Comp.2的贡献率为23.22%，第三主成分Comp.3的贡献率为11.59%，前三个主成分的方差贡献率累计为86.6%，大于80%或85%，因此选择此三个主成分即可较典型的表达该地区的农业生态经济的发展状况，从主成分示意图中亦可以看出Comp.3为图形由陡峭变为平坦的转折点，另外从主成分散点图中可以直观看出各原始变量对主成分的贡献情况。写出主成分与原变量的线性关系式为：
F1=0.342*X1-0.368*X2-0.375*X4-0.355*X5+0.312*X6+0.559*X7+0.113*X8+0.233*X9
0.614*X2+0.155*X4-0.761*X5-0.110*X6
F3=-0.446*X1+0.206*X6+0.467*X7-0.203*X8-0.692*X9
根据主成分与原始变量之间的线性关系式，给出各成分的实际解释意义：
a. 第一主成分F1是多个变量的线性组合，且各变量的系数大小相当，与X1、X6、X7呈现较强的正相关，与X2、X4、X5、X9呈现较强的负相关，这几个变量综合反映了生态经济结构状况，因此我们认为第一主成分F1是生态经济结构的总体代表；
b. 第二主成分F2与X2、X4正相关，与X5呈现较强的负相关，而X2、X4、X5分别为人均占有资源情况，因此认为第二主成分F2代表了人均资源量；
c. 第三主成分F3与X2、X8呈现较强的负相关，与X6、X7、X9呈现正相关， 这几个变量体现的是土地面积的使用效率，因此认为第三个主成分代表了土地使用情况；
d. 显然，用三个主成分F1、F2、F3代替原来9个变量（X1~X9）来描述该地区农业生态经济系统，使得问题更简化，增强了对变量的概括性解释。
10.1 健康实例1：在埃及一个叫做卡拉马的村庄，测量了161个儿童的身高，年龄和身高数据见下表，试分析年龄和身高之间的关系，并预测年龄为30个月时的身高。
年龄（月）
平均身高（厘米）
height=c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8,82.8,83.5)
plot(age,height)
lm.reg=lm(formula=height~age)
summary(lm.reg)
abline(lsfit(age,height))
lm(formula = height ~ age)
Residuals:
-0.248 -0.014
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
& 2e-16 ***
29.66 4.43e-11 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.256 on 10 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic:
880 on 1 and 10 DF,
p-value: 4.428e-11
& age=age[-3];height=height[-3]
& lm.reg2=lm(formula=height~age)
& summary(lm.reg2)
lm(formula = height ~ age)
Residuals:
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
& 2e-16 ***
35.64 5.33e-11 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2073 on 9 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic:
1270 on 1 and 9 DF,
p-value: 5.327e-11
& age.pre=data.frame(age=30)
& h.pre=predict(lm.reg2,age.pre,interval=&prediction&,level=0.95)
lwr
upr
1 84.02034 83.46839 84.57228
（1）首先，为了直观显示年龄和身高之间的关系，先画出一张散点图（上图左），从图中看见年龄和身高基本在一条直线附近，可以认为两者具有线性关系；然后建立回归模型，回归分析结果为：y=64.x；之后，进行F检验，得到P=4.428e-11远小于0.05，说明线性模型的系数非常显著。
（2）接下来， 对回归方程的残差进行分析，画出残差图，发现第3个点的残差较大，有些异常；剔除异常值后，重新拟合回归方程，得到：y=64.554+0.6489x，F检验的P=5.327e-11变得更小了，说明取得第3个点之后，回归方程更显著了（有意义了）。
（3）根据上面回归方程，对年龄30个月时的身高进行预测估计，在0.95置信水平下，预测值为84.02034，预测区间为[83.28]。
10.2 健康实例2：根据人类生病率数据库（HMD,2007)中瑞典的人口信息，选取其中的男性数据，试对生病率和年份、年龄之间的关系，拟合模型并对未来进行预测。
Female_Exp
Female_dth
L_female_exp
L_male_exp
993.005341
108.9732528
120.0060706
29.9985862
11.0847134
& setwd(&d:/R311/data_23714&)
& mort=read.csv('swedish_mortality.csv',header=T)
& mort=na.omit(mort)
& mort=mort[mort$q_male&0,]
& mort=mort[mort$q_male&1,]
& attach(mort)
The following objects are masked from mort (pos = 3):
Age, Female_death, Female_Exp, L_female_exp, L_male_exp,
Male_death, Male_Exp, q_female, q_male, Year
The following objects are masked from mort (pos = 5):
Age, Female_death, Female_Exp, L_female_exp, L_male_exp,
Male_death, Male_Exp, q_female, q_male, Year& library(rgl)
& plot3d(Year,Age,q_male,col='grey',type='p',zlim=1)
& par(mfrow=c(1,2))
& plot(Age,log(q_male),col='blue',main=&年龄与生病率（对数）&)
& plot(Year,log(q_male),col='blue',main=&年份与生病率（对数）&)
& layout(1)
& plot(Age,L_male_exp,col='blue',main=&年龄与对数生存人数&)
& hist(Male_death,col='blue',freq=F,breaks=100) & lg.md=log(Male_death) & data1=lg.md[lg.md&6] & data2=lg.md[lg.md&6] & p1=length(data1)/(length(data1)+length(data2)) & library(MASS) & para1=fitdistr(data1,'normal')$estimate& para2=fitdistr(data2,'normal')$estimate& p=p1*dnorm(lg.md,para1[1],para1[2])+(1-p1)*dnorm(lg.md,para2[1],para2[2])& hist(log(Male_death),col='blue',freq=F,breaks=100,main=&双峰分布--经验值和拟合值&,ylim=c(0,0.33))& points(lg.md,p,col=2)& F.mix=p1*pnorm(lg.md,para1[1],para1[2])+(1-p1)*pnorm(lg.md,para2[1],para2[2])& ks.test(lg.md,F.mix)
Two-sample Kolmogorov-Smirnov test
lg.md and F.mix
D = 0.99, p-value & 2.2e-16
alternative hypothesis: two-sided
& ols=lm(Male_death~Age+Year+L_male_exp)& summary(ols)
lm(formula = Male_death ~ Age + Year + L_male_exp)
Residuals:
-485.0 -277.6
201.6 1727.7
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept) -2.66e+03
6.6e-07 ***
& 2e-16 ***
L_male_exp
& 2e-16 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 324 on 5722 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 3.04e+03 on 3 and 5722 DF,
p-value: &2e-16
& m.nb=glm.nb(Male_death~factor(Age)+factor(Year)+offset(L_male_exp))& summary(m.nb)& anova(m.nb,test='Chisq')
Analysis of Deviance Table
Model: Negative Binomial(114), link: log
Response: Male_death
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(&Chi)
factor(Age)
&2e-16 ***
factor(Year)
&2e-16 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.
& par(mfrow=c(2,2))& plot(m.nb,col='blue',main=&广义线性模型诊断图&)& location=c(1,102,307,) & mort=mort[-location,]& attach(mort)& model.final=glm.nb(Male_death~poly(Age,25)+poly(Year,4)+offset(L_male_exp)) & options(digits=2)& summary(model.final)& pre.final=predict(model.final) & layout(1) & plot(Male_death,exp(pre.final),col='blue',xlab='观测值',ylab='拟合值',main='正交多项式改进模型')& abline(0,1,col='red')& q_pre=exp(pre.final)/Male_Exp & plot(Age,log(q_male),col='blue',pch='.',main='对数生病率') & points(Age,log(q_pre),pch='.',col=2)& legend(70,-7,legend=c('观察值','拟合值'),lty=1,col=c(1,2))
a. 整理数据。在瑞典年生病率数据集中，共包含6000多条记录，部分记录有缺失，使用na.omit()函数剔除不完整记录。对于生病率而言，取值范围为0&q&=1，不在此范围的点均视为异常值，也要剔除。
b. 绘制散点图。散点图是查看待研究变量之间关系的最直观方式。三维、二维散点图如上图所示。从中可以发现以下规律：
随着年龄增加，男性人口生病率总体趋势上升，但在0~20岁阶段，死亡率由新生儿时期的高值逐渐降低，随后升高，并具有一定波动性，呈现“浴盆曲线”。人口老龄化近年来越发严重，生病率（对数）随着时间的推移不断降低。生存人数随着年龄升高逐渐减少。
c. 人口数已知，生病率建模转化为生病人数建模。从直方图中可以看出来，死亡人数是一个右偏分布，右尾拖得很长，很难找到一个合适的分布去拟合它。于是对变量取对数，再进行拟合。
d. 取对数后，生病人数呈现“双峰”分布，对于双峰分布，一般采用分段拟合的方法。假设数据由两个均值不同、方差不同的正态分布混合而成，选择图形中以横轴数值6为分段界线。分段拟合后，其拟合曲线与直方图基本一致。
e. 用KS检验来判断双峰分布的拟合效果是否有效，计算得p-value
& 2.2e-16远小于0.05，应该拒绝原假设，说明分段拟合的双峰分布模型及参数还无法确切描述生病人数模型，需要探索其它模型。
f. 考虑到获得的数据集中包含多个变量，可以从变量关系的角度构建模型，而研究多个变量关系，最简单的方式是普通多元线性回归模型。用lm()拟合后得Year的P值0.8，大于0.05，说明普通多元线性回归模型并不理想。
g. 考虑广义线性模型。将每一个年龄、每一个年份都作为一个因子水平，整个模型中共包含164个参数（109个年龄+54个年份+1原始），选择负二项广义线性模型glm.nb()。其拟合结果中，几乎所有因子（只有几个例外）水平都显著（即P值小于0.05），初步认为拟合效果还行，下面进行诊断并改进模型。
h.下面用方差分析anova()，检验各个因子的显著性。
计算结果中，两个因子Age、Year都高度显著(P&2e-16远小于0.05），说明它们对生病人数的估计占有重要作用。最后，绘制残差的散点图进行广义线性模型的拟合效果诊断。
i. 残差诊断图共有四组，第一张是残差的散点图，数据点基本均匀地分布在横轴y=0两侧，诊断出两个异常点——样本307和5626；第二张是QQ图，图中的点大部分都集中在y=x直线附近，同样诊断出两个异常点；第三张是位置-尺度图，样本的残传值越小，散点分布越靠下，而位置较高的几个点则为异常点；第四张是曲式距离图，曲式距离表示每一个点对回归模型的影响力，各样本点应比较均匀，曲式距离偏大的点就可以认为是异常点。于是，剔除四图中标识出的异常点。
j. 剔除异常值后，还需要控制模型中参数个数，因为参数越多，估计结果误差越大。减少模型参数的一个快捷办法是在回归方程中对年龄Age和年份Year分别引入正交多项式，用poly()实现，并用AIC准则作为判断标准。经尝试发现，poly(Age,25)+poly(Year,4)时，AIC最小（AIC=53272），从而得知30个参数的最优模型（25+4+1）。
k. 此时再用负二项广义线性模型拟合，发现尽管残差deviance有所增加，但AIC统计量明显改善（变小）。
选择此时的模型为最终模型。
l. 重新查看拟合效果。用predict()计算拟合值，画出图形，可以看出散点集中于直线y=x，样本观测值与拟合值比较一致，模型拟合（预测）效果较好。
m. 以上建模找到了生病人数与年份、年龄之间的关系，而生病人数除以生存人数即为生病率，而与生病率最密切的变量是年龄（而非年份），所以需要找到生病率随年龄变化的模型，然后对年龄增大后生病率进行预测，见最后一张图形。
11.1 财政实例：财政收入是一个国家经济系统最重要的部分，通过分析财政一般收入中的各分项，能够及时发现现行税制和政策是否适应经济发展情况、产业结构是否合理等问题。财政收入会受到各种不同因素的影响，为深入研究，以财政收入y（亿元）为因变量，自变量选取如下：第一产业国内生产总值x1（亿元）、第二产业国内生产总值x2（亿元）、第三产业国内生产总值x3（亿元）、人口数x4（万人）、社会消费品零售总额x5（亿元）、受灾面积x6（万公顷），试建立计量经济学模型。
& revenue=read.table(&d:/R311/xxf_R/xxf_data/revenue.txt&,header=T)
& lm.reg=lm(y~x1+x2+x3+x4+x5+x6,data=revenue)
& summary(lm.reg)
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = revenue)
Residuals:
-295.71 -173.52
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
18.829 8.12e-11 ***
-1.171e-01
15.067 1.31e-09 ***
-5.152e-01
2.930e-02 -17.585 1.91e-10 ***
-1.104e-01
0.00206 **
-1.864e-02
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 234.8 on 13 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 2.294e+04 on 6 and 13 DF,
p-value: & 2.2e-16
3）中间结果：
利用多元线性模型来进行回归计算，拟合优度R^2=0.9999，说明模型的拟合效果很好，但在多元情况下因为自变量个数较多，拟合优度会自然较高，所以还要看检验的结果；回归方程的F检验P&2.2e-16远远小于0.05，所以线性模型是合适的，只是回归系数x1、x2的对应P值0.107不小0.05，即不显著，回归系统x6对应的P值0.0显著性水平下不显著，在0.1显著性水平下才显著，于是决定从模型中剔除自变量x1、x2，但留下其它自变量，也留下x6自变量（也可以考虑剔除X6，本案例选择保留），重新进行拟合分析。
4）再次分析：
& lm.reg1=update(lm.reg,.~.-x1-x2)
& summary(lm.reg1)
lm(formula = y ~ x3 + x4 + x5 + x6, data = revenue)
Residuals:
-325.62 -147.54
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
27.020 3.89e-14 ***
& 2e-16 ***
-5.438e-01
1.981e-02 -27.445 3.09e-14 ***
-1.392e-01
-7.256 2.80e-06 ***
-1.803e-02
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 233.6 on 15 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 3.476e+04 on 4 and 15 DF,
p-value: & 2.2e-16
从再次分析的结果中，可知x3、x4、x5对应的P值均远小于0.05，x6的P值也小于0.1（虽然不小于0.05），可以认为自变量系数均较为显著，写出拟合后的回归方程如下：Y=4(X3)-0...01803(X6)，即财政收入受第三产业国内生产总值、人口数量、社会消费品零售总额、受灾面积的影响，且具有多元线性关系。
12. 消费实例：我国2003年各地区农村居民家庭平均每人主要食品消费量如表所示，为了在市场分析中根据调查变量建立指标体系，请对数据进行综合分析。
6蛋类及制品
& food=read.table(&d:/R311/data_23714/food.txt&,header=T)
& food=food[,-1]
& library(labdsv)
& food.pca=pca(food,dim=4,cor=TRUE)
& summary(food.pca)
Importance of components:
Standard deviation
1.1 1.21
Proportion of Variance 0.7 0.05
Cumulative Proportion
0.9 0.33
& loadings.pca(food.pca)
＿ -0.550
＿ -0.242
＿ -0.320 -0.616
＿ 0.186 -0.697 -0.470
＿ -0.168
0.509 -0.142
＿ -0.263
0.408 -0.270 ＿
0.119 -0.112 -0.367
＿ 0.201
& op=par(mfrow=c(1,2))
& varplot.pca(food.pca)
从方差贡献率Proportion of Variance和直方图中可知，前四个主成分贡献率累计：32.52%+24.67%+14.93%+8.90%=81.02%，大于80%，可用前四个主成分来代替原始的9个变量；第一个主成分PC1与X4~X8相关性较高（系数较大），说明第一主成分主要反映猪牛羊肉、家禽、水产品等方面，可以归为肉制品类；第二个主成分PC2与X1、X8、X9相关性较高，分别对应粮食、食糠、酒，可以归为粮食类；第三个主成分PC3与X2、X3相关性较高，反映蔬菜、食油类；第四个主成分PC4与各原始变量的相关性分配比较平均，可以认为PC4反映的是所有食品的综合情况。据此，将原始变量分类，建立市场分析中调查研究时的指标体系。
13.1 营销案例1：在市场营销领域，营销研究往往结合市场调查问卷进行，以了解消费者需求动态，分析产品综合竞争力，为经营管理决策提供可靠的理论依据。某公司想要了解消费者购买牙膏时更追求什么目标，于是通过商场拦访对30人进行了访谈，用7级“克里特”量表（1表示非常不同意，7表示非常同意）询问他们对以下陈述的认同程度：
V1——购买预防蛀牙的牙膏是重要的；
V2——我喜欢使牙齿亮泽的牙膏；
V3——牙膏应当保护牙龈；
V4——我喜欢使口气清新的牙膏；
V5——预防坏牙不是牙膏提供的一项重要功效；
V6——购买牙膏时最重要的考虑是如何创造富有魅力的牙齿。
& yagao=read.table(&d:/R311/data_23714/yagao.txt&,header=T)
& data=yagao[,-1]
& factanal(data,factors=2)
factanal(x = data, factors = 2)
Uniquenesses:
0.104 0.440 0.141 0.384 0.238 0.294
Factor1 Factor2
0.945 ＿
＿ 0.747
＿ 0.779
V5 -0.868 ＿
＿ 0.838
Factor1 Factor2
SS loadings
Proportion Var
Cumulative Var
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 4.86 on 4 degrees of freedom.
The p-value is 0.302
本案例中有6个变量V1~V6，首先设定提取因子（主成分）数为2，用factanal()函数进行主成分分析后，观察到方差贡献率累计值为73.3%。通常要求方差贡献率累计达到80%以上，但本案例中原始变量个数较少（只有6个），用少数几个因子很难提取出大量的原始信息，因此70%以上的贡献率也可以接受（但若太低，就需要增加所设定的因子个数）。另外，发现主成分（因子）F1、F2对应的变量系数区分明显，故认为将因子个数选择为2是有效且合理的。2个主成分（因子）与原变量之间的线性关系如下：
0.945*V1+0.917
*V3-0.868*V5
0.747*V2-0.135*V3+0.779 *V4+0.838*V6
解释实际意义：
a. F1和变量V1、V3、V5相关性很强，分别对应“预防蛀牙”、“保护牙龈”、“预防坏牙”，因此F1可以看做是保健利益因子；
b. F2和变量V2、V4、V6相关性很强，分别对应“牙齿亮泽”、“口气清新”、“富有魅力”，因此F2可以看做是社交利益因子；
c. 通过拦访、问卷调查、数据分析，可以得出结论，消费者想从牙膏中得到的利益有两大类：保健利益和社交利益，供经营管理者决策参考。
13.2 营销案例2：某企业对218名顾客进行了收入水平和品牌选择关系的调查研究，收入水平分为高、中、低3个等级，品牌有A~F6个，数据见下表，试分析收入水平和品牌选择之间的关联性。
& brand=data.frame(low=c(2,49,4,4,15,1),medium=c(7,7,5,49,2,7),high=c(16,3,23,5,5,14))
& rownames(brand)=c(&A&,&B&,&C&,&D&,&E&,&F&)
& library(ca)
& options(digits=3)
& brand.ca=ca(brand)
& brand.ca
Principal inertias (eigenvalues):
Percentage 60.75%
0.266 0.101
1.029 0.738
0.282 0.055
-0. -0.581 -0.850 0.988 -0.8296
1.368 -1.403 0.281
low medium
& names(brand.ca)
&rownames&
[6] &rowinertia& &rowcoord&
&colnames&
[11] &coldist&
&colinertia& &colcoord&
& brand.ca$rowcoord
[1,] -0.727
1.399 -0.200
[3,] -0.581
[4,] -0.850 -1.403
[6,] -0.830
& plot(brand.ca,main=&收入水平与品牌选择之间的关联性示意图&)
从关联分析（对应分析）示意图中的横坐标距离关系可以看出，低收入人群倾向于选择品牌B和E，中收入水平人群倾向于选择品牌D，而高收入水平倾向于选择品牌A、C和F，此分析结果可以帮助企业进行市场定位和营销决策。
13.3 营销案例3：企业对市场上的产品进行分类，从而更准确的制定营销策略。某饮料企业收集了市场上16种饮料的热量、咖啡因、钠含量和价格4种指标数据，见下表，试进行市场细分。
& drink=read.table(&d:/R311/data_23714/drink.txt&,header=T)
& drink=drink[,-1]
& d=dist(drink)
& hc1=hclust(d,method=&ward.D&)
& hc2=hclust(d,method=&single&)
& hc3=hclust(d,method=&complete&)
& hc4=hclust(d,method=&median&)
& opar=par(mfrow=c(2,2))
& plot(hc1,hang=-1);plot(hc2,hang=-1);plot(hc3,hang=-1);plot(hc4,hang=-1)
& par=(opar)
& cutree(hc1,4)
[1] 1 2 2 2 3 3 3 2 3 1 2 4 2 4 3 3
& cutree(hc2,4)
[1] 1 2 2 2 2 2 3 2 2 1 2 4 2 4 2 2
& cutree(hc3,4)
[1] 1 2 3 2 4 3 4 2 3 1 2 2 2 2 4 4
& cutree(hc4,4)
[1] 1 2 2 2 3 3 3 2 3 1 2 4 2 4 3 3
采用4种不同的方法（离差平方和法、最短距离法、最长距离法、中间距离法）将16种饮料进行分类，因方法不同，所得到的结果也略有出入，其中最短距离法与其余三种差别最大。观察示意图确定要分为4类较合适（或者分成5类也可以），然后用序号1、2、3、4分别表示4个类别，实现市场细分。
14.1 经济案例1：随着电子信息技术的发展，以电话和快递等为代表的邮电业在近几年发展迅速，越来越多的人选择网上或电话购物，以快递的形式发送或接收货物。邮电业的发展与第三产业的发展密不可分，必须适应我国当前的经济结构，为了找出邮电业和经济发展的深层次关系，根据《中国统计年鉴》的国民经济数据，来衡量邮电业发展情况和经济发展情况，数据见下表。
信函（亿件）
快递（万件）
移动电话年末用户（万户）
固定电话年末用户（万户）
第一产业总产值（万亿元）
工业总产值（万亿元）
建筑业总产值（万亿元）
第三产业总产值（万亿元）
& post=read.table(&d:/R311/data_23714/post.txt&,header=T)
& post=post[,-1]
& post=scale(post) #对数据做正态离差标准化
& ca=cancor(post[,1:4],post[,5:8])
& options(digits=4)
[1] 0.2 0.7
0.00 -0.012
X2 -0.02 -0.667
X4 -0.13 -1.337
Y1 -0.03 -0.4
Y3 -0.05 -9.0
Y4 -0.26 -0.7
5.871e-18 -6.192e-17
7.686e-17 -9.714e-17 -1.183e-16
从相关系数$cor可以看出，前两队数据U1和V1、U2和V2之间具有高度正相关关系（r1=0.9984，r2=0.9512），说明邮电业越发达，经济形势越好。因此可以根据这两队变量来分析邮电业和经济发展的相关关系。前两对标准化的典型变量的线性组合(U1,V1)、(U2,V2)是：
U1=0.027*X1-0.054*X2-0.133*X3-0.119*X4
V1=-0.086*Y1+0.065*Y2-0.059*Y3-0.209*Y4
U2=-0.18*X1-0.246*X2+1.659*X3-1.513*X4
V2=-0.117*Y1+1...261*Y4
a. 在第一对典型相关变量（U1,V1）中，U1是邮电业各指标的线性组合，其中X3（移动电话年末用户）和X4（固定电话年末用户）比其他变量的系数更大，说明移动电话和固定电话是邮电业的主要指标，在邮电业中占主导地位；而快递比函件的系数要大，说明随着经济的发展，快递对邮电业的贡献在增加，函件在减少；V1是经济发展各指标的线性组合，其中Y4（第三产业增加值）的系数最大，说明第三产业是各产业中与邮电业相关联的主要指标，而其他变量系数基本相当，说明它们与邮电业的关联性都差不多；
b. 在第二对典型相关变量（U2,V2）中，U2是邮电业各指标的线性组合，其中仍是X3（移动电话年末用户）和X4（固定电话年末用户）的系数最大；固定电话和移动电话用户量的符号相反，说明它们之间有一定的相互抑制作用；V2是经济发展各指标的线性组合，其中Y2（工业）、Y3（建筑业）、Y4（第三产业）的系数较大，根据它们的符号判断，Y2（工业）、Y3（建筑业）和X3（移动电话年末用户）是同方向发展的，Y4（第三产业）和X4（固定电话年末用户）是同方向发展的。
14.2 经济案例2：自金融危机后，世界经济从2009年开始逐渐复苏，我国的经济水平也延续了一直以来的高速发展，然而伴随而来的是物价飞涨，利用我国2005年1月至2012年11月共119期宏观数据（CPI消费者价格指数、M2广义货币、NEER人民币BIS名义有效汇率、PPI工业品出厂价格指数、IAV工业增加值指数、API农产品价格指数），试分析CPI上涨的原因及传导机制，建立模型对CPI未来走向做出短期预测。
& data=read.table(&d:/R311/data_23714/cpi_data.txt&,header=T)
& dat=ts(data,frequency=12,start=c(2005,1))
& dat=log(dat)
& plot.ts(dat,main=&CPI时间序列&,col='blue')
& dat_dif=diff(dat)
& plot.ts(dat_dif^2,col='blue') & library(tseries) & cpi=ts(data$CPI,frequency=12,start=c(2005,1))
& cpi_dif=diff(log(cpi))
& pp.test(cpi_dif)
Phillips-Perron Unit Root Test
Dickey-Fuller Z(alpha) = -113, Truncation lag parameter = 4, p-value =
alternative hypothesis: stationary
& install.packages(&vars&)
& library(vars)
& options(digits=2)
& result=VAR(dat_dif,p=2)
& result$varresult$CPI
& summary(result)$varresult$CPI
lm(formula = y ~ -1 + ., data = datamat)
Residuals:
-0.....016608
Coefficients:
Estimate Std. Error t value Pr(&|t|)
NEER.l1 -0.099707
0.00096 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0062 on 103 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 2.97 on 12 and 103 DF,
p-value: 0.00137
& dat_mod=dat_dif[,c(-2,-4)]
& result2=VAR(dat_mod,p=2)
& summary(result2)$varresult$CPI
lm(formula = y ~ -1 + ., data = datamat)
Residuals:
-0.....017203
Coefficients:
Estimate Std. Error t value Pr(&|t|)
NEER.l1 -0.079505
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0061 on 107 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 4.14 on 8 and 107 DF,
p-value: 0.000247
& plot(result2)
等待页面改变的确认...
等待页面改变的确认...
等待页面改变的确认...
& pre=predict(result2,n.ahead=12)
& CPI.pre=pre$fcst$CPI
lower upper
4.2e-03 -0. 0.012
[2,] -3.4e-04 -0. 0.013
4.8e-04 -0. 0.014
[4,] -9.4e-05 -0. 0.014
6.5e-04 -0. 0.014
2.8e-04 -0. 0.014
1.4e-04 -0. 0.014
9.2e-05 -0. 0.014
1.2e-04 -0. 0.014
1.2e-04 -0. 0.014
1.1e-04 -0. 0.014
1.0e-04 -0. 0.014
& options(digits=5)
& LCPI.pre=log(cpi[119])+cumsum(CPI.pre[,-1])
& CPI.pre=exp(LCPI.pre)
[1] 101.204
98.990 100.376 101.781 103.204 104.453 105.856 107.306
[28] 108.791 110.298 111.829 113.380 114.953 116.549 118.166 119.806 121.468
& plot(CPI.pre,col='blue',main=&对未来12期CPI数值进行预测&)
a. 向量自回归模型VAR是预测工具之一，VAR模型要求参加建模的变量满足平稳性。在对数据进行处理之前，对原始数据取对数，主要有以下好处：取对数可以将高级乘法关系变为线性关系；使得数量级变小，可以消除一定的序列中的异方差性、共线性、非平稳性等；减少数据处理过程中可能产生的误差。
b. 取对数后，观察数据的平稳性，发现6个变量几乎都不平稳，于是将变量逐一进行一阶差分，然后进行平稳性检验（PP检验），这里以6个变量中CPI为例，其PP检验结果p-value =
0.01，小于0.05,说明较显著，即一阶差分后平稳，同理可检验其它5个变量，结果为均一阶差分后平稳。
c. 以CPI为例进行VAR建模。从VAR建模结果的CPI相关Coeficients中可以看出各参数估计值及其显著性，其中M2和PPI对CPI十分不显著，即货币供应量、工业品出厂价格指数对CPI的解释性不强，因此M2和PPI不能作为CPI的解释变量。
d. 剔除M2和PPI之后，重新VAR建模。可以发现，剔除M2和PPI后，剩余参数更加显著。根据参数估计结果，写出VAR模型的第一个方程：
LogCPI=0.52*LogCPI(t-1)-0.07951*LogNEER(t-1)+0.01023*LogIAV(t-1)+0.12603*LogAPI(t-1)+0.22404*LogCPI(t-2)+0.07385*LogNEER(t-2)+0.00759*LogIAV(t-2)-0.02189*LogAPI(t-2)
由于t-2期的参数并不十分显著，这里可以主要分析t-1期系数所反映的因素影响程度。
首先对CPI影响最大的是API农产品价格指数，这与我国食品通胀严重的情况相符，农产品价格的增长直接带动食品价格增加，而我国通货膨胀的60%来自食品。
此外，IAV工业增加值也会带动CPI的同时变动，这与实际情况也是一致的。
人民币的实际汇率系数为负，原因有可能是增加的货币供给并未出现在流通领域，而是广泛地用于房地产投资、金融衍生品投资等领域，因此并未造成物价上涨，反而因为这种投资挤压了消费需求，使得物价下跌，正好印证了目前我国房地产投资过热现象。
e. 根据VAR模型拟合曲线和残差图可知，VAR模型大致可以描述经济指数时间序列的整体走势，但短期的拟合效果要好于长期；而残差值围绕0值上下波动，没有明显的变化趋势；从残差自相关图中可知，CPI、NEER、IAV、API的前5阶系数基本都落入置信界限内部（只有API的第3阶超出）。
f. 利用VAR模型进行预测。
第一列fcst是对未来12期的预测值，lower和upper是预测区间上下界；
g. 由于前文VAR模型是先对数再差分的数据结果，所以需要把预测值进行差分和对数的逆运算，获得最终预测结果。
15.1 企业案例1：在诊断中小企业经营状态时，经常选择4个经济指标用于判断企业处于破产状态还是正常运行状态，其中X1总负债率（现金收益/总负债）、X2收益性指标（纯收入/总财产）、X3短期支付能力（流动资产/流动负债）、X4生产效率性指标（流动资产/纯销售额）。已知17个破产企业（记为1类）和21个正常运行企业（记为2类）的数据资料，现对其它8个未知状态的企业进行判断。
已知状态企业
收益性指标
短期支付能力
生产效率性指标
待判状态企业
& B=read.table(&d:/R311/data_23714/bankruptcy.txt&,header=T)
& mu=colMeans(B)
& Sx=cov(B)
& distance=mahalanobis(B,mu,Sx)
& options(digits=3)
& distance
4.438 10.528
0.578 30.306
2.955 12.755
& library(WMDB)
& G=c(rep(1,17),rep(2,21))
& G=as.factor(G)
& wmd(B,G)
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
blong 1 1 1 1 1 1 1 1 2
27 28 29 30 31 32 33 34 35 36 37 38
[1] &num of wrong judgement&
9 11 12 14 25 30 32
[1] &samples divided to&
[1] 2 2 2 2 1 1 1
[1] &samples actually belongs to&
[1] 1 1 1 1 2 2 2
Levels: 1 2
[1] &percent of right judgement&
& newdata=data.frame(X1=c(0.04,-0.06,0.07,-0.13,0.15,0.16,0.29,0.54),
X2=c(0.01,-0.06,-0.01,-0.14,0.06,0.05,0.06,0.11),
X3=c(1.5,1.37,1.37,1.42,2.23,2.31,1.84,2.33),
X4=c(0.71,0.4,0.34,0.44,0.56,0.2,0.38,0.48))
& wmd(B,G,TstX=newdata)
1 2 3 4 5 6 7 8
blong 1 1 1 1 2 2 1 1
根据马氏距离的判别分析可知，8个待判企业中，除企业5和企业6处于正常运行状态外，其余6个企业均处于破产状态。
15.2 企业案例2： 某企业整理出了近年来各季度的销售数据，如下表，试从中分析销售额的动态趋势，以便更好的指导生产。
& data=c(362,385,432,341,382,409,498,387,473,513,582,474,
+ 544,582,681,557,628,707,773,592,627,725,854,661)
& sales=ts(data,frequency=4,start=c(2004,1))
& plot.ts(sales)
& components=decompose(sales)
& options(digits=3)
& components$seasonal
73.4 -64.8
73.4 -64.8
73.4 -64.8
73.4 -64.8
73.4 -64.8
73.4 -64.8
& plot(components)
-------------------------------------------------------------------------
& data=c(362,385,432,341,382,409,498,387,473,513,582,474,544,582,681,557,628,707,773,592,627,725,854,661)
& sales=ts(data,frequency=4,start=c(2004,1))
& sales.pre=HoltWinters(sales)
& sales.pre
Holt-Winters exponential smoothing with trend and additive seasonal component.
HoltWinters(x = sales)
Smoothing parameters:
alpha: 0.3612309
beta : 0.1084233
Coefficients:
s1 -30.83908
s3 147.48133
s4 -66.32662
& plot(sales.pre,main=&用HoltWinters模型对企业销售额数据进行拟合&)
& library(forecast)
& sales.pre2=forecast.HoltWinters(sales.pre,h=4)
& plot(sales.pre2,main=&用HoltWinters模型对未来1年（4个季度）销售额数据进行预测&)
& acf(sales.pre2$residuals,lag.max=8)
& Box.test(sales.pre2$residuals,lag=8,type=&Ljung-Box&)
Box-Ljung test
sales.pre2$residuals
X-squared = 9.4184, df = 8, p-value = 0.3082
a. 从“Time-sales”图形可以看出，本企业的销售额有着明显的趋势性（向上发展）和季节性，需要将趋势性和季节性分解开来，如图“Time-observed-trend-seasonal-random”，可以看出每年前三个季度销售额上升，而第四个季度下降。
b. 利用HoltWinters模型对销售额数据进行拟合，从图中可以看出红线和黑线基本一致；再利用HoltWinters模型对未来1年（4个季度）的销售额进行预测，曲线末端的浅灰色阴影表示90%置信水平的预测区间，深灰色阴影表示80%置信水平的预测区间，中间蓝点表示预测值。
c. 利用自相关图对HoltWinters模型拟合（预测）的销售额的残差进行检验，发现各阶（0阶除外）系数均落入自相关置信界限，所以可初步判断残差无明显相关性；再利用Ljung-Box模型检验残差，检验结果p-value = 0.3082，大于0.05显著性水平，没有足够理由拒绝原假设，说明残差不存在相关性，即认为HoltWinters模型对本案例数据的拟合（预测）效果很好。.
16. 天气案例：已知年100年间伦敦每年的降雨量（单位：英尺），数据来源：。试分析此时间序列的规律特征，并对未来10年降雨量进行预测。
& data=scan(&/tsdldata/hurst/precip1.dat&,skip=1)
Read 100 items
& rain=ts(data,start=1813)
& plot.ts(rain,main=&年100年间伦敦每年的降雨量数据（单位：英尺）&)& rain.pre=HoltWinters(rain,beta=FALSE,gamma=FALSE)
& rain.pre
Holt-Winters exponential smoothing without trend and without seasonal component.
HoltWinters(x = rain, beta = FALSE, gamma = FALSE)
Smoothing parameters:
beta : FALSE
gamma: FALSE
Coefficients:
a 24.67819
& plot(rain.pre)
& rain.pre$SSE
& library(forecast)
& rain.pre2=forecast.HoltWinters(rain.pre,h=10)
& plot.forecast(rain.pre2,col=2)
& acf(rain.pre2$residuals,lag.max=20)
a. 从降雨量曲线图可知，降雨量在25英尺左右的水平上下波动，并未有明显的趋势性和季节性，可以认为其是大致平稳的，而随机波动在整个时间序列期间也没有明显的增加或减小，所以可以用加法模型来描述此数据的时间序列特征，而数据的预测应选用简单指数平滑法。
b. 经HoltWinter模型拟合后，原始降雨量数据的波动量较大曲线变成一条较为平缓的曲线，此拟合曲线不包含随机波动，主要反映平稳时间序列的恒定水平，用此简单指数平滑法做拟合的样本内拟合误差（即残差）的平方和SSE为。
c. 继续用HoltWinter模型对未来十年降雨量进行预测，
图中浅灰色阴影区为置信水平95%的预测区间，深灰色阴影区为置信水平80%的预测区间，中间粗线条为预测的降水量水平。
d. 判断时间序列模型拟合效果的重要指标是残差，残差属于样本内的预测误差，若残差属于白噪声序列，则说明“残差彼此独立、不存在显著的自相关、预测效果较好”。绘制伦敦降雨量延迟20阶的残差自相关图，图中两条虚线表示置信界限，属于自相关系数的上下限。若残差自相关系数超过边界（0阶处除外），则表明存在相关关系。而本降雨量案例中，拟合（预测）残差的所有阶（0阶除外）自相关系数都落入了置信界限内，一般可以断定残差序列不存在自相关，拟合（预测）效果较好。
17.1 投资案例1：证券投资组合是指投资者对各种证券商品进行一定的选择，从而形成相对固定的若干个投资品种，以达到在一定的约束下，实现投资收益最大化的目标。投资组合最优化问题，就是寻找何种证券组合才是最有效的投资组合。选择上证50概念股票中的前30支股票（按股票代码顺序）在日0时前220个交易日的日收益率数据，使用不考虑现金红利的日个股回报率，试分析最优化投资组合。
& setwd(&d:/R311/data_23714&)
& data=read.csv(&Dalyr.csv&,header=TRUE)
& code=unique(data$Stkcd)
& n=length(code)
& Dalyr=matrix(0,220,n)
& for(i in 1:n)
+ Dalyr[,i]=data$Dretnd[which(data$Stkcd==code[i])]
& r.cov=cov(Dalyr)
& r.mean=apply(Dalyr,2,mean)
& options(digits=2)
[1] -0.059 -0.016 -0.015 -0.056
[10] -0.056
0.046 -0.028
0.043 -0.00189
[19] -0.038
0.037 -0.010 -0.00113
& plot(r.mean,main=&30支股票期望收益率&,pch=8,col=2)
& ###朴素方法：均值-方差模型###
& t11=Sys.time()
& num=万次循环，画出随机点
& rp=numeric(num); #定义变量，存放投资组合的均值
& sigmap=numeric(num) #定义变量，存放投资组合的标准差
& for(i in 1:num)
+ x1=runif(30)
+ x=x1/sum(x1)
+ rp[i]=sum(r.mean*x)
+ sigmap[i]=t(x)%*%r.cov%*%x
& plot(sigmap,rp,pch='.',main=&随机打点法的有效边界&)
& ###计算有效边缘
& rp=round(rp,4)
& rp.uni=unique(rp)
& nn=length(rp.uni)
& sigma.min=numeric(nn)
& for(i in 1:nn)
+ sigma.min[i]=min(sigmap[which(rp==rp.uni[i])])
& rp.sort=sort(rp.uni)
& order=order(rp.uni)
& lines(rp.sort~sigma.min[order],col=&red&)
& t12=Sys.time()
& time.used1=t12-t11
& time.used1
Time difference of 6.6 secs
-------------------------------------------------------------------------------
& ###模拟退火算法（simulated annealing）###
& t21=Sys.time()
& temp0=1000 #确定初始“温度”
& temp=temp0
& x0=c(1,runif(29))
& x=x0/sum(x0)
& M=1000 #惩罚因子
& R=0.0025 #目标收益率
& f1=function(x){t(x)%*%r.cov%*%x+M*max(0,R-sum(r.mean*x))}
& f2=function(x){(R-sum(r.mean*x))/sqrt(t(x)%*%r.cov%*%x)}
& f=as.numeric(f)
#记录退温次数
& a=data.frame(x)
& e=0.00001
& while(temp&0.00001)
+ for(i in 1:L)
+ x1=x+runif(30)
+ x1=x1/sum(x1)
+ deltaf=f1(x1)-f1(x)
+ pr=exp(-deltaf/temp)
+ if(deltaf&0)
+ f[k+1]=f1(x1)
+ if(pr&runif(1))
+ f[k+1]=f1(x1)
+ temp=temp0/(k^3)
& #退火法选取的投资组合#
& r.sa=0;sigma.sa=0
& for(i in 1:length(a))
+ r.sa[i]=sum(r.mean*a[[i]])
+ sigma.sa[i]=t(a[[i]])%*%r.cov%*%a[[i]]
& plot(sigma.sa,r.sa,main=&模拟退火法有效边界&)
& #有效边缘
& r.sa=round(r.sa,4)
& rsa.uni=unique(r.sa)
& nn=length(rsa.uni)
& sigma.min=numeric(nn)
& for(i in 1:nn)
+ sigma.min[i]=min(sigma.sa[which(r.sa==rsa.uni[i])])
& rsa.sort=sort(rsa.uni)
& order=order(rsa.uni)
& lines(rsa.sort~sigma.min[order],col='red', lwd=2)
& legend(0.e-04,legend=&正态假设法&)
& t22=Sys.time()
& time.used2=t22-t21
& time.used2
Time difference of 2.2 secs
a. 计算30支股票在220个交易日的平均收益率，绘制散点图；然后使用朴素（“笨”）方法画出股票投资组合的有效边界，即依据均值-方差模型的原理，随机模拟10万种投资组合的权重，将所有投资组合的点画在均值方差图中，选取每个期望收益下风险最小的点连成有效边界曲线，共耗时6.6s。但这种方法过于朴素（思路简单，当计算系统能力强大时可以采用。），对于现实中远多于30支股票的较大计算量，需要找到更好的组合最优化算法，例如：模拟退火算法，可以节省时间，但思路复杂。
采用退火法计算并画出投资组合的有效边界，从而节省时间，但由于退火法的循环次数过少，使得模拟的有效载荷边界不完整，且投资组合点显得过于分散。
17.2 投资案例2：从网站获取多支股票的交易数据，然后构建资产组合的有效前沿。
1）数据：来自网络，实时获取，见分析代码部分。
& library(quantmod)
& getSymbols(c('IBM','SPY','YHOO'))
& IBM_ret=dailyReturn(IBM)
& SPY_ret=dailyReturn(SPY)
& YHOO_ret=dailyReturn(YHOO)
& data=merge(IBM_ret,SPY_ret,YHOO_ret)
& library(timeSeries)
& data_ts=as.timeSeries(data)
利用getSymbols()从网络实时获取股票数据，然后可以用“17.1 投资案例1”中的方法获得投资组合的最优前沿，也可以用R2.9版本中的fPortfolio程序包来快速画出投资组合最优前沿，具体如下：
library(fPortfolio)
Frontier=portfolioFrontier(data_ts)
plot(Frontier)
但是，fPortfolio程序包不能应用于R2.9之后的版本之中，若你的R版本大于2.9，那么你需要重新安装R2.9版本，然后再利用fPortfolio程序包。
17.3 投资案例3：先确定投资组合，再确定组合的有效边界！前面两个投资案例提到了股票投资组合的最优化问题，即针对股票组合寻找有效边界，那么市场上股票有几千支，评价指标很多（侧重不同），如何选择投资组合呢？途径就是：把股票指标综合到一起对股票进行分类。本案例选取钢铁行业的20支股票及7个指标作为分类的原始数据和依据。
总资产（亿元）
主营收入（亿元）
净利润增长率
每股净资产（元）
净资产收益率
主营业务收入增长率
每股资本公积金（元）
2）分析： & stock=read.table(&d:/R311/data_23714/stock.txt&,header=T)
& rownames(stock)=stock[,1]
& KM=kmeans(stock[,2:8],4)
K-means clustering with 4 clusters of sizes 5, 1, 2, 12
Cluster means:
107.0 598.292 4..028
2 9. 4..320
471.0 682.040 4..228 129.800
136.5 248.960 3..540
Clustering vector:
大钢不锈 安阳钢铁 鲁银投资 南钢股份 武钢股份 菜钢股份 柳钢股份 凌钢股份
华凌股份 济南钢铁 唐钢股份 杭钢股份 安泰股份 承德钒钛 韶钢松山 本钢板材
八一钢铁 宝钢股份
鹏博士 广钢股份
Within cluster sum of squares by cluster:
(between_SS / total_SS =
Available components:
[1] &cluster&
&withinss&
[5] &tot.withinss& &betweenss&
[9] &ifault&
& sort(KM$cluster)
安阳钢铁 鲁银投资 南钢股份 菜钢股份 柳钢股份 宝钢股份 大钢不锈 武钢股份
凌钢股份 华凌股份 济南钢铁 唐钢股份 杭钢股份 安泰股份 承德钒钛 韶钢松山
本钢板材 八一钢铁
鹏博士 广钢股份
------------------------------------------------------------------------------------------------------
& subset=subset(stock,select=2:8)
& d=dist(subset)
& hc=hclust(d,method=&ward.D&)
& #plclust(hc)
& plot(hc)
& rect.hclust(hc,k=4,border=&red&)
首先利用kmesns算法对钢铁股票信息数据进行分类，分成了1、2、3、4四个类别；然后利用hclust()层次聚类分析进行分类，分类结果见以上2图，二者的结果有时会有一点不完全一致。第一类属于传统大规模钢铁联合企业，综合竞争力很强，这样的股票，无论涨或跌，幅度都不会很大；第三类属于总体能力处于均衡的位置，有一定投资效果，但不明显；第四类属于潜力型股票，未来前景一片大好。然后，就可以从股票分类结果中，选择股票组合，然后再进行组合的有效边界计算，然后选择最优情况进行投资。
参考文献：
1. 《数据分析：R语言实战》，李诗羽，张飞，王正林 编著。
-----------------------------------------------------------------------------------------------------------------------------------
为您业务中的“疑难杂症”提供针对性解决方案，是我们永恒的追求！
Mobile: 187
（若网页内容，存在不合适之处，亦恳请您告诉我们，谢谢！）
