##########Duplan sztochasztikus betamodell. IID exponencialisok a csucsokba. Exponencialis parameterenek ##########becslese #betamodelles gráf generálása betagraf=function(x){ n=length(x) betagraf=sapply(c(1:n), function(i){rbinom(n,1,x[i]*x/(1+x[i]*x))}) betagraf[lower.tri(betagraf)]=t(betagraf)[lower.tri(betagraf)] diag(betagraf)=0 betagraf } #adott gráf valószínûsége logaritmusának meghatározása vsz=function(x,graf){ n=length(x) seged=sapply(c(1:n), function(i){x[i]*x/(1+x[i]*x)}) seged[lower.tri(seged)]=1 diag(seged)=1 seged2=1-seged seged2[lower.tri(seged2)]=1 diag(seged2)=1 vsz=sum(log(seged^graf*seged2^(1-graf))) } ##gráf fokszáma egyszerûen a rowSums ##Metropolis-Hastings logkülönbség. x1 és x0 a csúcsokra írt számok logkul=function(x1,x0,graf,landa){ A1=x1%*%t(x1) diag(A1)=0 A0=x0%*%t(x0) diag(A0)=0 logkul=0.5*sum(log(1+A0)-log(1+A1))+sum(rowSums(graf)*(log(x1)-log(x0)))+landa*sum(x0-x1) } EMexp=function(graf,n1,n2,n3,landa0,x0,szigma){ m=length(x0) landa=landa0 a=1 for (j in 1:n3){ if (abs(a)>10^(-3)){ for (i in 1:n1){ x1=pmax(10^(-3),x0+rnorm(m,sd = szigma)) if (logkul(x1,x0,graf,landa)>log(runif(1))){ x0=x1 } } y=0 for (i in 1:n2){ x1=pmax(10^(-3),x0+rnorm(m,sd = szigma)) if (logkul(x1,x0,graf,landa)>log(runif(1))){ x0=x1 } y=y+sum(x0) } a=m*n2/y-landa landa=m*n2/y } } landa } print(date()) eredlanda_1=sapply(c(1:2), function(i){EMexp(betagraf(rexp(1100)),100,1000,30,0.5,rep(1,1100),0.1)}) print(date()) ##########Veletlen matrix meret 800x800tol 1000x1000ig 100asaaval lepkedve ###########Kell parallel csomag eredlanda_1=matrix(0,100,3) #eredmenyek matrixa library(parallel) print(date()) for (j in 0:2){ print(j) print(date()) no_cores=detectCores()-1 cl=makeCluster(no_cores) clusterExport(cl,"betagraf") clusterExport(cl,"vsz") clusterExport(cl,"logkul") clusterExport(cl,"EMexp") clusterExport(cl,"eredlanda_1") clusterExport(cl,"j") eredlanda_1[,(j+1)]=parSapply(cl,c(1:100),function(i){EMexp(betagraf(rexp((800+j*100))),100,1000,30,0.5,rep(1,(800+j*100)),0.1)}) print(date()) stopCluster(cl) } save(eredlanda_1,file="eredlanda_1.RData") print(date())