symulacje

Zastosowanie algorytmu Gibbsa

N=1000
lambda=matrix(0,N,10)
czas=c()
beta=c()
# wprowadzenie danych
dane=c(5,1,5,14,3,19,1,1,4,22)
czas=c(94.32,15.72,62.88,125.76,5.24,31.44,1.05,1.05,2.10,10.48)
alfa=1.8
beta=c(); beta[1]=0
delta=1
gamma=0.01

for (i in 2:N)
{
 for (j in 1:10)
 lambda[i,j]=rgamma(1,sh=dane[j]+alfa,ra=czas[j]+beta[i-1])
 beta[i]=rgamma(1,sh=gamma+10*alfa,ra=delta+sum(lambda[i,]))
}
# podział okna wykresu
# 1 wiersz i 3 kolumny
par(mfrow=c(1,3))
# wyświetlenie wyników
hist(lambda[,1])
hist(lambda[,2])
hist(beta)