jetzt kann MCS laufen

sim/theo.R

@@ -5,8 +5,9 @@ M              <- 100  # Anz. Messpunkte pro SZ
 N              <- 5    # Anz. der SZ
 p.mean         <- 0.0  # mbar
 m.mean         <- -1.8e-3 # mbar/s
+t.between      <- 5.0
 t.rand.sd.rel  <- 1e-3
-p.rand.sd.rel  <- 1e-3
+p.rand.sd.rel  <- 1e-8
 p.mean.sz      <- rep(p.mean, N) # Druckmittelwert der SZ
 m.sz           <- rep(m.mean, N) # Ausgangssteigung der SZ
 idx            <- seq(1, M)
@@ -27,18 +28,21 @@ p.rand.sd      <- abs(p.fill * t.rand.sd.rel)
 
 p.new  <- matrix(NA, ncol = M , nrow = N)
 p.rand <- matrix(NA, ncol = M , nrow = N)
+p.rand <- matrix(NA, ncol = M , nrow = N)
 t.new  <- matrix(NA, ncol = M , nrow = N)
 t.rand <- matrix(NA, ncol = M , nrow = N)
 
-m.rand <- rep(NA, N)
-c.rand <- rep(NA, N)
+m.rand       <- rep(NA, N)
+c.rand       <- rep(NA, N)
+p.rand.mean  <- rep(NA, N)
+t.rand.mean  <- rep(NA, N)
 
 for(i in 1:N){
 
   if(i == 1){
     t.new[i,]  <- t.d[i] * idx
   }else{
-    t.new[i,]  <- t.new[i-1, M] + t.d[i] * idx
+    t.new[i,]  <- t.new[i-1, M] + t.between + t.d[i] * idx
   }
 
   p.new[i,]  <- seq(p.upper[i],
@@ -48,20 +52,25 @@ for(i in 1:N){
   p.rand[i,] <- p.new[i,] + rnorm(M, 0, p.rand.sd)
   t.rand[i,] <- t.new[i,] + rnorm(M, 0, t.rand.sd)
 
-  temp.lm     <- lm( p.rand[i,] ~ t.rand[i,])
-  c.rand[i]   <- coefficients(temp.lm)[1]
-  m.rand[i]   <- coefficients(temp.lm)[2]
-
+  temp.lm        <- lm( p.rand[i,] ~ t.rand[i,])
+  c.rand[i]      <- coefficients(temp.lm)[1]
+  m.rand[i]      <- coefficients(temp.lm)[2]
+  p.rand.mean[i] <- mean(p.rand[i,])
+  t.rand.mean[i] <- mean(t.rand[i,])
 }
 
+c.i      <- p.rand.mean - m.rand * t.rand.mean
+t.rand.0 <- (mean(p.rand.mean) - c.i)/ m.rand
+
+Delta.t <- diff(t.rand.0)
+
 pt.df <- data.frame(p.new  = as.vector(p.new),
                     t.new  = as.vector(t.new),
                     p.rand = as.vector(p.rand),
                     t.rand = as.vector(t.rand))
-
-
 ggplot(pt.df) +
   geom_point(aes(x = t.new, y = p.new), color='black') +
   geom_point(aes(x = t.rand, y = p.rand), color='red') +
   geom_abline(aes(intercept = c.rand, slope = m.rand), color='red') +
+  geom_vline(xintercept = t.rand.0) +
   ggtitle("dp vs. dt")