希望能向大神们在R语言编程方面多多请教~~


@任坤
@谢益辉
  1. R语言在国内的应用水平还不高,大多还处在一般应用级别!可参见我2018年初所写的《2017年R语言发展报告(国内)》

2017年R语言发展报告(国内)?

mp.weixin.qq.com图标

2. 顶部的R玩家(高手)国内还有不少,汤银才老师、谢益辉、张丹、陈堰平、肖楠、刘思喆等很多高手。

可以多关注《R语言中文社区》,关注各领域大牛的技术实践分享


@Song

欢迎关注我们啊,我们写过很多关于R语言的文章,当然,存货更多,有疑问也必定知无不言。一起交流啊~~


关于昆仑数据:昆仑数据是工业大数据领域的领军企业,创始团队来自国内顶尖信息科技与工业企业以及顶尖研究机构,蝉联「中国大数据企业50强」,受邀参与制订《中国制造2025》工业大数据技术路线图,发起成立并主导运营工业大数据制造业创新中心,致力于用大数据和人工智慧技术。已服务新能源、石油天然气、电子制造、工程机械、环保、动力装备、生物制药等领域。更多内容请关注微信公众号 ( id:k2datas ) 了解。

编辑于 2018-03-23继续浏览内容知乎发现更大的世界打开Chrome继续statoorstatoor玩数如厨艺,预测如算命@肖凯@TomHall@统计之都
@肖凯@TomHall@统计之都
必须@自己一下,哈哈。关于R和并行计算的,可以看下ParallelR网站。http://www.parallelr.com/blog
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)


有很多搞生物搞计算的都是大神,比如@Charles Liang。

@HopeR


推荐阅读:
相关文章