筋肉神経接合シュミレーション
前回のMIKUで話題になったシュミレーション
自分にも出来るかRでやってみた。
# ベクトルpのindexを筋肉の番号、数値を神経の番号とする
# １回の単位時間が経つと全ての神経が違う筋肉支配を一つしようと試みる
# 新たな支配が上手くいくかは、その時における神経の支配筋肉数に比例するとする
N<- 100
p<- sample(N, replace = F)
p
# まず１回目の試行を考えてみる
p1<- sample(N, replace = F)
for( i in 1 : N ){
a<-p[ i ]
b<-p1[ i ]
num.a <- length( p[ p=a ] )
num.b <- length( p[ p=b ] )
prob.i <- c( num.a , num.b )/(num.a + num.b )
new.i <- sample(c( a ,b ) , prob = prob.i )
p[ i ] <- new.i [1]
}
for( j in 1:100){
p1<- sample(N, replace = F)
for( i in 1 : N ){
a<-p[ i ]
b<-p1[ i ]
num.a <- length( p[ p=a ] )
num.b <- length( p[ p=b ] )
prob.i <- c( num.a , num.b )/(num.a + num.b )
new.i <- sample(c( a ,b ), prob = prob.i )
p[ i ] <- new.i [1]
}
}
# あれ？もっと階段状になる予定…
# しかも、そんなに入れ替わってないぞ？
# 原因；new.i <- sampleのところ？ pのindexの取り出し方が間違い
for( j in 1:100){
p1<- sample(N, replace = F)
for( i in 1 : N ){
a<-p[ i ]
b<-p1[ i ]
num.a <- length( p[ p==a ] )
num.b <- length( p[ p==b ] )
prob.i <- c( num.a , num.b )/(num.a + num.b )
new.i <- sample(c( a ,b ), (num.a + num.b ) , prob = prob.i , replace=T)
p[ i ] <- new.i [1]
}
}
plot(sort(p))
#ぽいのが出来た！
# あとは分かりやすくplotしよう
# １の数は何個？
length( p[ p==1 ] )
num.p <- rep(0, N )
#これを全体に適応ー＞新しいベクトルに個数を代入
for( i in 1: N){
num.i <- length( p[ p==i ] )
num.p[i] <- num.i
}
# こんな感じ？N <-1000でやってみる
N <- 1000
for( j in 1:100){
p1<- sample(N, replace = F)
for( i in 1 : N ){
a<-p[ i ]
b<-p1[ i ]
num.a <- length( p[ p==a ] )
num.b <- length( p[ p==b ] )
prob.i <- c( num.a , num.b )/(num.a + num.b )
new.i <- sample(c( a ,b ), (num.a + num.b ) , prob = prob.i , replace=T)
p[ i ] <- new.i [1]
}
}
num.p <- rep(0, N )
for( i in 1: N){
num.i <- length( p[ p==i ] )
num.p[i] <- num.i
}
plot(sort(num.p ))
＃空間条件を入れだすともうちょっとややこしくなりそうです
