Computerintensive Methoden - Übung 1

From StatWiki

pdf(rpdf)
y<-c(2.5, 3.2, 5.0, 2.2, 3.2, 6.0)
m<-c(0,1,1,0,0,1)

daten<-data.frame(cbind(y,m))


gibbs.qtl<-function(daten,it) {

			# it...Anzahl der Iterationen
			# daten...erste Spalte ist der Phänotyp, zweite Spalte der Marker

r=0.1
sigma=1


mu0<-vector(length=it+1)
mu1<-vector(length=it+1)

mu0[1]=1
mu1[1]=3

for (i in 2:(it+1)) {

q<-vector(length=dim(daten)[1])


	for (j in 1:dim(daten)[1]) {

		if(daten[j,2]==0) {
			p<-dnorm(daten[j,1],mean=mu0[i],sd=sigma)*(1-r)/(dnorm(daten[j,1],mean=mu0[i],sd=sigma)*(1-r)+dnorm(daten[j,1],mean=mu1[i],sd=sigma)*r)
		}
		if(daten[j,2]==1) {
			p<-dnorm(daten[j,1],mean=mu0[i],sd=sigma)*r/(dnorm(daten[j,1],mean=mu0[i],sd=sigma)*r+dnorm(daten[j,1],mean=mu1[i],sd=sigma)*(1-r))
		}

		q[j]<-rbinom(1,1,prob=p)
	
	}	
	
	y0<-mean(daten[,1][q==0])
	y1<-mean(daten[,1][q==1])

	mu0[i]<-rnorm(1,mean=y0,sd=sigma/sum(q==0))
	mu1[i]<-rnorm(1,mean=y1,sd=sigma/sum(q==1))
}


return(cbind(mu0,mu1))
}


erg<-gibbs.qtl(daten,1000)

plot(erg[,1],erg[,2])

#

Meines:

d=data.frame(y=c(2.5,3.2,5,2.2,3.2,6), m=c(0,1,1,0,0,1))

m0 = 1 m1 = 3 z2 = 1 r = 0.1

res = t(replicate(10, { run = apply(d, 1, function(row) { y = row[1] m = row[2]

# rekombination d0 = dnorm(y, m0, z2) d1 = dnorm(y, m1, z2) p = d0 * (1-r)^m / (d0 * (1-r)^m + d1 * r^m)

q = rbinom(1, 1, p)

if(q == 0) c(0, y, 0, 1) else c(y, 0, 1, 0) })

ym = rowMeans(run[1:2,]) nm = rowSums(run[3:4,])

m0 <<- rnorm(1, ym[1], z2/nm[1]) m1 <<- rnorm(1, ym[2], z2/nm[2])

c(m0, m1) }))

plot(V1~V2, data=res)