2.1 对数正态分布不完全k阶矩公式

LN<-function(mu,sigma,L,U,k){
Lk<-(log(L)-mu)/sigma-k*sigma
Uk<-(log(U)-mu)/sigma-k*sigma
exp(k*mu+k^2*sigma^2/2)*(pnorm(Uk)-pnorm(Lk))
}

verify the pdf

LN(8.5,0.8,0,Inf,0)

exam 2.1

EX=LN(8.5,0.8,0,Inf,1)
EX1=0.75*EX
EX2=LN(8.5,0.8,0,25000,1)+25000*LN(8.5,0.8,25000,Inf,0)

exam 2.2

VARX=LN(8.5,0.8,0,Inf,2)-LN(8.5,0.8,0,Inf,1)^2
VARX1=0.75^2*VARX
EX2_SQ=LN(8.5,0.8,0,25000,2)+25000^2*LN(8.5,0.8,25000,Inf,0)
VARX2=EX2_SQ-EX2^2
VARX2

2.2 正态分布不完全k阶矩公式

2.2.0 正态分布不完全0阶矩公式:CDF

pnorm(Inf)-pnorm(0)

2.2.1 正态分布不完全1阶矩公式

norm<-function(mu=0,sigma=1,L,U){
L1<-(L-mu)/sigma
U1<-(U-mu)/sigma
mu*{pnorm(U1)-pnorm(L1)}-sigma*{dnorm(U1)-dnorm(L1)}
}
norm(,,0,1e10)
norm(,,-Inf,Inf)

2.2.2 正态分布不完全2阶矩公式

norm2<-function(mu=0,sigma=1,L,U){
L1<-(L-mu)/sigma
U1<-(U-mu)/sigma
(mu^2+sigma^2)*{pnorm(U1)-pnorm(L1)}-
sigma*{(2*mu+sigma*U1)*dnorm(U1)-(2*mu+sigma*L1)*dnorm(L1)}
}
norm2(,,0,Inf)

/*按照变量date的值来拆分数据集
如date=19960131生成一个数据集，date=19960228生成一个
现在有很多个，到20080731了
*/
/*方法一
思路是先把distinct date的值变为宏变量，然后用data步输出。
反复读数据，会很慢
*/
data a;
do date=1991 to 2011;
output;
end;
run;
data null;
set a end=last;
call symputx('date'||left(n),date);
if last then call symput ('n',n);
run;
%macro split;
%do i=1 %to &n;
data b&i;
set a;
where date=&&date&i;
run;
%end;
%mend;
%split;

/方法二/
data date;
date=input("19960131",yymmdd8.);
output;
date=input("19960228",yymmdd8.);
output;
date=input("19960309",yymmdd8.);
output;
run;

data null;
set date;
call execute ("data data_"||strip(put(date,yymmddn8.))||";set a;if date="||date||";run;");
run;

/方法三/

data a;
format target_day $8.;
input target_day $;
cards;
20120730
20120731
20120810
20120817
;run;

%macro segmentation();
proc sort data = a out = work.t1(keep = target_day) nodupkey;
by target_day;
run;
data null;
set work.t1 nobs=n;
call symputx(compress("set_"||n), put(input(target_day, yymmdd8.), yymmddn8.));
if n = 1 then call symputx('n',n);
run;
data %do ii = 1 %to &n.;
data_&&set_&ii..
%end;;
set a;
select (target_day);
%do ii = 1 %to &n.;
when ("&&set_&ii..") output data_&&set_&ii..;
%end;
otherwise;
end;
run;
proc sql;
drop table work.t1;
quit;
%mend;
%segmentation();

3 罚广义线性模型 51

install.packages("lars")
install.packages("ridge")

(1) OLS回归

library(ridge)
data(GenCont)
model.ols <- lm(Phenotypes ~ ., data = as.data.frame(GenCont))
summary(model.ols)

提取回归系数

coef.ols <- coef(model.ols)

查看不等于0的回归系数

coef.ols[coef.ols!=0]  

(2) ridge regression（岭回归）

model.rid <- linearRidge(Phenotypes ~ ., seq(0, 10, length = 100), data = as.data.frame(GenCont))

查看结果

summary(model.rid)

提取系数

coef返回的是标准化系数; 直接输入model.rid返回的是原始数据系数
coef.rid <- coef(model.rid)

model.rid2 <- linearRidge(Phenotypes ~ ., data = as.data.frame(GenCont))
summary(model.rid2)
coef.rid <- coef(model.rid2)

查看不等于0的系数

coef.rid[coef.rid!=0]                           

(3) 做 lasso regression

模型设定

library(lars)
dt<-as.matrix(GenCont)
X<-dt[,-1]
Y<-dt[,1]
model.lasso1 <- lars(X, Y)    
model.lasso2 <- lars(X, Y, type='lasso') 
model.lasso3 <- lars(X, Y, type='for')  # Can use abbreviations

看变量选入顺序

model.lasso3
summary(model.lasso3)

Cp的含义：衡量多重共线性，其取值越小越好，这里取到第3步使得Cp值最小

画个图

plot(model.lasso3)
set.seed(12345)

做交叉验证

CV.lasso <- cv.lars(X, Y, K=10)   

(best <- CV.lasso$index[which.min(CV.lasso$cv)])

选择最好的效果

(coef.lasso <- coef.lars(model.lasso3, mode='fraction', s=best))

命名

names(coef.lasso) <- colnames(dt)[-1]

查看结果

coef.lasso[coef.lasso!=0]