setwd("F:/Users/Miklos/Documents/Kutatas/TG/betamodel/Cikk") #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:100), function(i){EMexp(betagraf(rexp(30)),100,1000,30,0.5,rep(1,30),0.1)}) print(date()) eredlanda_05=sapply(c(1:100), function(i){EMexp(betagraf(rexp(30, rate = 0.5)),100,1000,30,0.5,rep(1,30),0.1)}) print(date()) eredlanda_2=sapply(c(1:100), function(i){EMexp(betagraf(rexp(30, rate = 2)),100,1000,30,0.5,rep(1,30),0.1)}) print(date()) print(date()) eredlanda_landa_1_meret_15=sapply(c(1:100), function(i){EMexp(betagraf(rexp(15)),100,1000,30,0.5,rep(1,15),0.1)}) print(date()) eredlanda_landa_1_meret_60=sapply(c(1:100), function(i){EMexp(betagraf(rexp(60)),100,1000,30,0.5,rep(1,60),0.1)}) print(date()) eredlanda_landa_1_meret_100=sapply(c(1:100), function(i){EMexp(betagraf(rexp(100)),100,1000,30,0.5,rep(1,100),0.1)}) print(date()) kidob=sapply(c(1:10), function(i){EMexp(betagraf(rexp(15)),100,1000,30,0.5,rep(1,15),0.1)}) print(date())