数理神経生物学 J.S.グリフィス著をもとにプログラムを作ってみました
library(dplyr)
lambda=0.3
exp_f=function(x){
z<-x*lambda*exp(-lambda*x)
return(z)
}
#平均発火時間
ave_t=integrate(exp_f,0,Inf)
#時刻０から時刻ｔの間に細胞がｘ回発火する確率
pois=function(x,t){
z<-(exp(-lambda*t)*(lambda*t)^x)/factorial(x)
return(z)
}
pois(3,2)
#指数分布を正規乱数として、系列相関係数の値により、独立かどうか調べる(kendall,1948)
samples=rexp(1000,rate=lambda)
n=length(samples);
r_data=data.frame(p=1:(length(samples)-1),r=0)
for(j in 1:(nrow(r_data)-1)){ r_data$r[r_data$p==p]
p=r_data$p[j]
samples_data=data.frame(num=1:length(samples),sample=samples)
A00=sum(samples_data$sample[1:(n-p)]^2)/(n-p)
A00_matrix=array(0,dim=c(n-p,n-p))
for(l in 1:(n-p)){
for(k in 1:(n-p)){
A00_matrix[l,k]=samples[l]*samples[k]
}
}
A00=A00-sum(A00_matrix)/((n-p)^2)
A0p=sum(samples_data$sample[1:(n-p)]*samples_data$sample[(p+1):(n)])/(n-p)
A0p_matrix=array(0,dim=c(n-p,n-p))
for(l in 1:(n-p)){
for(k in 1:(n-p)){
A0p_matrix[l,k]=samples[l]*samples[k+p]
}
}
A0p=A0p-sum(A0p_matrix)/((n-p)^2)
App=sum(samples_data$sample[(p+1):(n)]^2)/(n-p)
App_matrix=array(0,dim=c(n-p,n-p))
for(l in 1:(n-p)){
for(k in 1:(n-p)){
App_matrix[l,k]=samples[l+p]*samples[k+p]
}
}
App=App-sum(App_matrix)/((n-p)^2)
r_data$r[r_data$p==p]=A0p*(A00*App)^(-1/2)
}
plot(r_data$p,r_data$r,xlab="p",ylab="r")
#ｔ検定により、H0:average=0を検定(生物統計学入門;p.74)
r=r_data$r
u=sd(r);
alpha=0.05
if(abs(mean(r))<qt(1-alpha/2,nrow(r_data)-1)*(u/sqrt(nrow(r_data)))){
print("Hypothesis H0 is accepted")
}
#random walk(確率p,qに従いそれぞれ細胞にインパルスが到達したり、しなかったりする)
p=0.6;q=1-p
#到達点ｍに行くには、正の方向をa、負の方向をbとしたとき、a-b=m,a+b=t;合計遷移時刻にならなければいけない
arrived_dis=function(t,m){
if(floor((t+m)/2)==(t+m)/2){
z<-(factorial(t)/(factorial((1/2)*(t+m))*factorial((1/2)*(t-m))))*(p^((1/2)*(t+m)))*(q^((1/2)*(t-m)))
z[is.nan(z)==T]=0
}else{
z<-0
}
return(z)
}
#random walk の例
m=10;
times=9
N=100
for(j in 1:times){
random=sample(c(-1,1),N,prob=c(p,q),replace = T)
random=cumsum(random)
if(j!=times){
if(length(random[random>m])>0){
plot(c(1:N),random,xlab="t",ylab="a+b",xlim=c(1,N),ylim=c(0,m),col=j)
par(new=T)
}
}else{
plot(c(1:N),random,xlab="t",ylab="a+b",xlim=c(1,N),ylim=c(0,m),col=j)
}
}
#到達確率のプロット
library(pforeach)
prob_data=pforeach(j=1:100,.c=rbind)({
data=data.frame(m=j,t=1:100,prob=0)
for(l in 1:nrow(data)){
data$prob[l]=arrived_dis(data$t[l],data$m[l])
}
data
})
prob_data=prob_data %>% dplyr::filter(prob!=0)
library(rgl)
plot3d(prob_data$m,prob_data$t,prob_data$prob,type="s",col=2,xlab="m",ylab="t",zlab="prob")
rm(t,m)
#許容された経路の確率(m<sita)
path_mean=function(t,m,sita){
if(floor((t+m)/2)==(t+m)/2){
z<-((factorial(t)/(factorial((1/2)*(t+m))*factorial((1/2)*(t-m))))-(factorial(t)/(factorial((1/2)*(t-m)+sita)*factorial((1/2)*(t+m)-sita))))*(p^((1/2)*(t+m)))*(q^((1/2)*(t-m)))
}else{
z<-0
}
return(z)
}
path_mean(10,4,5)
#(16)式より関数を定義
I=function(t,sita){
if(sita!=0 & t!=0){
z<-path_mean((t-1),(sita-1),sita)*p
}else{
z<-1
}
return(z)
}
#(18)式より、関数を定義
I_re=function(t,sita){
z<-(sita/t)*(factorial(t)/(factorial((1/2)*(t+sita))*factorial((1/2)*(t-sita))))*(p^((1/2)*(t+sita)))*(q^((1/2)*(t-sita)))
return(z)
}
I_re(10,4);I(10,4)
#Sを定義
I_data=pforeach(j = 1:50,.c=rbind)({
data=data.frame(sita=j,t=c((j+1):150),I=0)
data
})
for(j in 1:nrow(I_data)){
I_data$I[j]=I_re(I_data$t[j],I_data$sita[j])
}
plot3d(I_data$sita,I_data$t,I_data$I,xlab="sita",ylab="t",zlab="I",col=2,type="s")
S=I_data %>% group_by(sita) %>% summarise(S=sum(I))
#sitaごとに計算
for(j in 1:50){
plot(I_data$t[I_data$sita==j],I_data$I[I_data$sita==j],xlab="t",ylab="I",col=2,type="o",main=paste0("sita:",j))
}
con_I=function(t,sita){
z<-(sita/sqrt(8*pi*p*q))*(t^(-3/2))*exp(-((sita-(2*p-1)*t)^2)/(8*t*p*q))
return(z)
}
con_I_data=pforeach(j = 1:500,.c=rbind)({
data=data.frame(sita=j,t=c((j+1):1000),I=0)
data
})
for(j in 1:nrow(con_I_data)){
con_I_data$I[j]=con_I(con_I_data$t[j],con_I_data$sita[j])
}
for(j in 1:500){
plot(con_I_data$t[con_I_data$sita==j],con_I_data$I[con_I_data$sita==j],xlab="t",ylab="I",col=2,type="o",main=paste0("con_I_sita:",j))
}
sum_con_I=function(s,sita){
func=function(t){
y<-(sita/sqrt(8*pi*p*q))*(t^(-3/2))*exp(-((sita-(2*p-1)*t)^2)/(8*t*p*q))
return(y)
}
a=0.001;H=s;n=100000
Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));
for(j in 1:(n-2)){
coefficient=ifelse(j-floor(j/2)*2==0,4,2)
Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))
}
return(Y)
}
ave_con_I=function(s,sita){
func=function(t){
y<-(sita/sqrt(8*pi*p*q))*t*(t^(-3/2))*exp(-((sita-(2*p-1)*t)^2)/(8*t*p*q))
return(y)
}
a=0.001;H=s;n=100000
Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));
for(j in 1:(n-2)){
coefficient=ifelse(j-floor(j/2)*2==0,4,2)
Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))
}
return(Y)
}
#熱方程式の陰解法（クランクニコルソン）
#藪下信 計算物理２ p,24
#時刻0で状態１にいると仮定する
N=100
k=0.013
t=0.1
x=0.01
p=(k*t)/(x^2)
A=B=C=array(0,dim=c(N,N))
vec=c(1,rep(0,N-1))
for(j in 1:(N-1)){
A[j,(j+1)]=1
C[j,(j+1)]=-1
}
A=t(A)+A;diag(A)=rep(-2,N)
C=t(C)+C;diag(C)=rep(2-2/p,N)
diag(B)=rep(-2/p,N)
T_vec=vec
times=1000
T_matrix=array(0,dim=c(N,times))
T_matrix[,1]=vec
for(j in 1:(times-1)){
T_vec=solve(A+B)%*%C%*%T_vec
T_matrix[,(j+1)]=T_vec
}
plot(T_vec)
plot(T_matrix[,50])
plot(T_matrix[,100])