Computerintensive Methoden - Übung 1

From StatWiki
Jump to: navigation, search
> 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)
Warning message:
In rnorm(1, mean = y0, sd = sigma/sum(q == 0)) : NAs produced
>
> 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)