希望能向大神們在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


推薦閱讀:
相关文章