自身の手作りプログラム（数理神経生物学における熱拡散方程式）

数理神経生物学　J.S.グリフィス著をもとにプログラムを作ってみました

```{r}

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)

  
}



```

```{r}
#熱方程式の陰解法（クランクニコルソン）

#藪下信　計算物理２ 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])

```