sum(pbino(22)*1:400) # 130.5237 for (i in 5:15) print(sum(pdarroch(n1,n2,i)*1:400)) # 279.6753 # 254.623 # 227.1532 # 199.8575 # 175.2202 # 154.4834 # 137.5962 # 123.8896 # 112.6351 # 103.2499 # 95.30768 sum(1:400*pcapture(2,n1+n2-m2,n1+n2)) # 168.4865 sum(1:400*pcapture(7,nplus,nc)) # 372.8943 g1=gibbs1(10000,7,nplus,nc,200) mean(g1$N) # 327.016 mean(g1$p) # 0.2268489 sq=seq(1,10000,by=20) par(mfrow=c(2,2)) plot(sq,g1$N[sq],ylab="N",xlab="",type="l",ylim=c(300,360)) plot(sq,g1$p[sq],ylab="p",xlab="",type="l",ylim=c(0.18,0.27)) hist(g1$N,nclass=100,xlim=c(300,360),prob=T,main="",xlab="N") hist(g1$p,nclass=100,xlim=c(0.18,0.27),prob=T,main="",xlab="p") dev.copy2eps(file="g1.eps") g2=gibbs2(10000,n1,c2,c3,200,10,5) mean(g2$N) # 57.0072 mean(g2$p) # 0.3985086 sq=seq(1,10000,by=20) par(mfrow=c(5,2),mar=c(2,4,1,1)) plot(sq,g2$p[sq],ylab="p",xlab="",type="l") hist(g2$p,nclass=100,prob=T,main="",ylab="") plot(sq,g2$q[sq],ylab="q",xlab="",type="l") hist(g2$q,nclass=100,prob=T,main="",ylab="") plot(sq,g2$N[sq],ylab="N",xlab="",type="l") hist(g2$N,nclass=100,prob=T,main="",ylab="") plot(sq,g2$r1[sq],ylab=expression(r[1]),xlab="",type="l") hist(g2$r1,nclass=100,prob=T,main="",ylab="") plot(sq,g2$r2[sq],ylab=expression(r[2]),xlab="",type="l") hist(g2$r2,nclass=100,prob=T,main="",ylab="") dev.copy2eps(file="g2.eps") plot(jitter(g2$r1,factor=1),jitter(g2$r2,factor=1),cex=0.5, xlab=expression(r[1]),ylab=expression(r[2])) dev.copy2eps(file="jit.eps") mean(seuil(0:11,n1,c2,c3,1,0.1)) # 0.1187102 ar1=ardipper(10000,n1,c2,c3,1,0.1) ar1=factor(ar1) barplot(summary(ar1)/10000,main="",xlab=expression(r[1])) dev.copy2eps(file="ar.eps") g3=gibbs3(10000,eurodip,eurodipabc) apply(g3$p,2,mean) # 0.2514232 0.2589981 0.2669728 apply(g3$phi,2,mean) # 0.9933285 0.9926140 0.9929496 mean(g3$psi[3,3,]) # 0.6089712 mean(g3$psi[1,2,]) # 0.3345478 sq=seq(20,10000,by=20) par(mfrow=c(3,2),mar=c(2,4,1,1)) plot(sq,g3$p[sq,1],ylab="p(1)",xlab="",type="l") hist(g3$p[20:10000,1],nclass=100,prob=T,main="",ylab="") plot(sq,g3$phi[sq,2],ylab=expression(phi(2)),xlab="",type="l") hist(g3$phi[20:10000,2],nclass=100,prob=T,main="",ylab="") plot(sq,g3$psi[3,3,sq],ylab=expression(psi(3,3)),xlab="",type="l") hist(g3$psi[3,3,20:10000],nclass=100,prob=T,main="",ylab="") dev.copy2eps(file="g3.eps")