フェノタイプの整理メモ１２ 確率的に考える７

• 病変は単発か多発か
  • 関連記事(こちら)
• 発病後、一定期間内の活動病変最大箇所数は、シミュレーションの条件により、複数個所をピークとするようになる

# 数えられるフェノタイプ

# 3ジェノタイプを比較する

NL<-1 # 座位数
TL<-120 # 一生は120年

# イベントがしょっちゅう起きて(均質化して)いる場合には
# Dxのタイミングにはばらつきが小さい

Fs<-0.01*2^(12) # イベントが起きる平均時間間隔(年)

# 少し逆行因子を強くする
Gs<-Fs*1.05 # イベントを修復する平均時間間隔(年)
ThresNums<-c(40)

Nlegion<-10 # 関節数とか

Niter<-100 # 人数

YearAfterDx<-10

NumActives<-rep(0,Niter)
MaxActives<-NumActives
for(fi in 1:length(Fs)){
	F<-Fs[fi]
	G<-Gs[fi]
	for(ti in 1:length(ThresNums)){
		ThresNum<-ThresNums[ti]
		for(ii in 1:Niter){
			NR<-ThresNum*1000
			Fevents<-matrix(rexp(Nlegion*NR,F),nrow=Nlegion)
			Gevents<-matrix(rexp(Nlegion*NR,G),nrow=Nlegion)
			FeventsCum<-t(apply(Fevents,1,cumsum))
			GeventsCum<-t(apply(Gevents,1,cumsum))
			
			FGtime<-cbind(FeventsCum,GeventsCum)
			FG11<-cbind(matrix(1,nrow=Nlegion,ncol=NR),matrix(-1,nrow=Nlegion,ncol=NR))
			zerotimesList<-c()
			zerotimes01List<-c()
			for(i in 1:Nlegion){
				ord<-order(FGtime[i,])
				tmp11<-FG11[i,ord]
				tmpcumsum<-cumsum(tmp11)
				truecumsum<-tmpcumsum
				
				negs<-which(tmpcumsum<0)
				if(length(negs)>0){
					negval<-tmpcumsum[negs]
					maxNeg<-min(negval)
					inits<-rep(0,abs(maxNeg))
					for(j in 1:length(inits)){
						inits[j]<-which(negval==-j)[1]
					}
					cancelled<-FG11[i,ord]
					cancelled[negs[inits]]<-0
					
					truecumsum<-cumsum(cancelled)
				}
				truecumsum01<-truecumsum-ThresNum
				zeros<-which(truecumsum01==0)
				if(length(zeros)>0){
					zerotimes<-FGtime[i,ord[zeros]]
					zerotimes01<-(-1)^(0:(length(zerotimes)-1))
					zerotimesList<-c(zerotimesList,zerotimes)
					zerotimes01List<-c(zerotimes01List,zerotimes01)
				}
				
				
#matplot(cbind(tmpcumsum,truecumsum),type="l")
			}
			ord2<-order(zerotimesList)
			#if(min(cumsum(zerotimes01List[ord2]))<0){
			plot(zerotimesList[ord2],cumsum(zerotimes01List[ord2]),type="s",xlim=c(0,TL))
			yearDx<-zerotimesList[ord2][1]
			NumTimePoints<-length(which(zerotimesList[ord2]<=yearDx+YearAfterDx))
			NumActives[ii]<-cumsum(zerotimes01List[ord2])[NumTimePoints]
			MaxActives[ii]<-max(cumsum(zerotimes01List[ord2])[1:NumTimePoints])
			#}
		}
	
	}

}

hist(NumActives,breaks=(-1):max(NumActives))
hist(MaxActives,breaks=(-1):max(MaxActives))
table(MaxActives)