用户名*
邮箱*
密码*
确认密码*
验证码* 点击图片更换验证码
找回密码
忘记密码了？输入你的注册邮箱，并点击重置，稍后，你将会收到一封密码重置邮件。
希望能向大神们在R语言编程方面多多请教~~
关于留学法国，实习工作，法国生活，欧洲旅游等问题，可发私信或加V：269511008
@任坤
@谢益辉
@Song
@肖凯@TomHall@统计之都
FinanceR – SegmentFault
最近一直在关注 一个叫FinanceR的专栏，作者更新速度很快，内容也是对的起业内良心呀。
链接在此，不要谢我！
最近在做一个模拟，出现了点问题，没有找出错误在哪儿，请教一下各位大神
##产生数据 generate.data=function(n,n0,c,Tpara) ##n0为非有效样本的个数 { x=runif(n,0,2) z=as.matrix(matrix(rbinom(2*n,1,0.5),nrow=n,ncol=2)) epsilon=rnorm(n,0,0.1) w=as.numeric((x+epsilon)>1) Tbeta=Tpara[1] Talpha=Tpara[2:3] lambda=(z%*%Talpha-1/sqrt(2))^2 h=exp(Tbeta*x+lambda) y=rexp(n,h) C=runif(n,0,c) T=apply(cbind(y,C),1,min) delta=as.numeric(y<=C) x1=sample(x,n0,replace=FALSE) ##非有效样本 eta=rep(1,n) for( i in 1:n0) { k=which(x==x1[i]) ##记录非有效样本的位置 eta[k]=0 } ##非有效样本的eta值为0,有效样本eta值为1 data0=data.frame(cbind(T,delta,x,w,eta,z)) data=data0[order(data0$T),] return(data) }
e.hat=function(beta,alpha,i) ## i表示从第i个个体开始累加 ,beta、alpha待估计 { a=eta*as.numeric(w==w[i])*exp(beta*x+z%*%alpha) b=eta*as.numeric(w==w[i]) num=sum(a[i:n]) ##分子 dem=sum(b[i:n]) ##分母 ehat=num/dem return(ehat) }
##Sfun实现一个向量从当前位置起后面元素累加,用到了函数cumsum,X为向量 Sfun=function(X) { d=sum(X)-cumsum(X) n=length(X) f=rep(0,n) g=d[-n] f[1]=sum(X) f[2:n]=g return(f) }
##似然函数 logLfun0=function(theta) { beta=theta[1] alpha=theta[2:3] e=rep(0,n) for(i in 1:n) { if(eta[i]==1) e[i]=0 else e[i]=e.hat(beta,alpha,i) } e=replace(e,is.na(e),1) r=exp(beta*x+z%*%alpha) c=eta*r+(1-eta)*e ##中间变量 L1=sum(delta*(log(c))) L2=sum(delta*log(Sfun(c))) ##Sfun为前面已经定义的函数 logL=L1-L2 ##对数似然函数 return(logL) }
##模拟M次 library(maxLik) main=function(M,n,n0,c,Tpara) ##求M次估计的平均值 { para=matrix(rep(0,3*M),M,3) for(i in 1:M ) { data=generate.data(n,n0,c,Tpara) T=data[,1] delta=data[,2] x=data[,3] w=data[,4] eta=data[,5] z=as.matrix(data[,6:7]) result=maxLik(logLfun0,start=c(runif(3,-1,1)),method=”NR”) para[i,]=round(result$estimate,2) } m=apply(para,2,mean)
s=apply(para,2,sd) #标准差
bias=Tpara-m #偏差
out=rbind(para,m,s,bias)
return(para)
}
M=20 ##模拟次数
n=100 ##每次模拟的样本量
n0=n*0.5 ##非有效样本的个数
c=2 ##删失比率大概为20%
Tpara=c(log(2),1,-0.5) ##真实参数，顺序为beta、alpha
如果我直接用命令：outcome=main(M,n,n0,c,Tpara)，R会报错：找不到对象eta.
如果我不把循环写在main函数里面而是直接用以下命令却没有问题，请假大神这是问什么？
para=matrix(rep(0,3*M),M,3)
for(i in 1:M )
{ data=generate.data(n,n0,c,Tpara)
T=data[,1]
delta=data[,2]
x=data[,3]
w=data[,4]
eta=data[,5]
z=as.matrix(data[,6:7])
result=maxLik(logLfun0,start=c(runif(3,-1,1)),method=”NR”)
para[i,]=round(result$estimate,2)
m=apply(para,2,mean)
s=apply(para,2,sd) ##M次估计的标准差，其中方差除以的是（n-1）
bias=Tpara-m
肖楠
必须@自己一下，哈哈。关于R和并行计算的，可以看下ParallelR网站。
@Yuxi Zhu大魔王
有很多搞生物搞计算的都是大神，比如@Charles Liang。
昵称*
E-Mail*
回复内容*
回复 ( 10 )
@任坤
@谢益辉
@Song
@肖凯@TomHall@统计之都
FinanceR – SegmentFault
最近一直在关注 一个叫FinanceR的专栏，作者更新速度很快，内容也是对的起业内良心呀。
链接在此，不要谢我！
最近在做一个模拟，出现了点问题，没有找出错误在哪儿，请教一下各位大神
##产生数据
generate.data=function(n,n0,c,Tpara) ##n0为非有效样本的个数
{ x=runif(n,0,2)
z=as.matrix(matrix(rbinom(2*n,1,0.5),nrow=n,ncol=2))
epsilon=rnorm(n,0,0.1)
w=as.numeric((x+epsilon)>1)
Tbeta=Tpara[1]
Talpha=Tpara[2:3]
lambda=(z%*%Talpha-1/sqrt(2))^2
h=exp(Tbeta*x+lambda)
y=rexp(n,h)
C=runif(n,0,c)
T=apply(cbind(y,C),1,min)
delta=as.numeric(y<=C)
x1=sample(x,n0,replace=FALSE) ##非有效样本
eta=rep(1,n)
for( i in 1:n0)
{ k=which(x==x1[i]) ##记录非有效样本的位置
eta[k]=0 } ##非有效样本的eta值为0,有效样本eta值为1
data0=data.frame(cbind(T,delta,x,w,eta,z))
data=data0[order(data0$T),]
return(data)
}
e.hat=function(beta,alpha,i) ## i表示从第i个个体开始累加 ,beta、alpha待估计
{ a=eta*as.numeric(w==w[i])*exp(beta*x+z%*%alpha)
b=eta*as.numeric(w==w[i])
num=sum(a[i:n]) ##分子
dem=sum(b[i:n]) ##分母
ehat=num/dem
return(ehat)
}
##Sfun实现一个向量从当前位置起后面元素累加,用到了函数cumsum,X为向量
Sfun=function(X)
{ d=sum(X)-cumsum(X)
n=length(X)
f=rep(0,n)
g=d[-n]
f[1]=sum(X)
f[2:n]=g
return(f)
}
##似然函数
logLfun0=function(theta)
{ beta=theta[1]
alpha=theta[2:3]
e=rep(0,n)
for(i in 1:n)
{ if(eta[i]==1)
e[i]=0
else
e[i]=e.hat(beta,alpha,i)
}
e=replace(e,is.na(e),1)
r=exp(beta*x+z%*%alpha)
c=eta*r+(1-eta)*e ##中间变量
L1=sum(delta*(log(c)))
L2=sum(delta*log(Sfun(c))) ##Sfun为前面已经定义的函数
logL=L1-L2 ##对数似然函数
return(logL)
}
##模拟M次
library(maxLik)
main=function(M,n,n0,c,Tpara) ##求M次估计的平均值
{ para=matrix(rep(0,3*M),M,3)
for(i in 1:M )
{ data=generate.data(n,n0,c,Tpara)
T=data[,1]
delta=data[,2]
x=data[,3]
w=data[,4]
eta=data[,5]
z=as.matrix(data[,6:7])
result=maxLik(logLfun0,start=c(runif(3,-1,1)),method=”NR”)
para[i,]=round(result$estimate,2)
}
m=apply(para,2,mean)
s=apply(para,2,sd) #标准差
bias=Tpara-m #偏差
out=rbind(para,m,s,bias)
return(para)
}
M=20 ##模拟次数
n=100 ##每次模拟的样本量
n0=n*0.5 ##非有效样本的个数
c=2 ##删失比率大概为20%
Tpara=c(log(2),1,-0.5) ##真实参数，顺序为beta、alpha
如果我直接用命令：outcome=main(M,n,n0,c,Tpara)，R会报错：找不到对象eta.
如果我不把循环写在main函数里面而是直接用以下命令却没有问题，请假大神这是问什么？
para=matrix(rep(0,3*M),M,3)
for(i in 1:M )
{ data=generate.data(n,n0,c,Tpara)
T=data[,1]
delta=data[,2]
x=data[,3]
w=data[,4]
eta=data[,5]
z=as.matrix(data[,6:7])
result=maxLik(logLfun0,start=c(runif(3,-1,1)),method=”NR”)
para[i,]=round(result$estimate,2)
}
m=apply(para,2,mean)
s=apply(para,2,sd) ##M次估计的标准差，其中方差除以的是（n-1）
bias=Tpara-m
out=rbind(para,m,s,bias)
肖楠
必须@自己一下，哈哈。关于R和并行计算的，可以看下ParallelR网站。
@Yuxi Zhu大魔王
有很多搞生物搞计算的都是大神，比如@Charles Liang。