"""对钙信号的动力学进行建模,AR模型。"""
import matplotlib.pyplot as plt
import numpy as np
if __name__ == '__main__':
length = 500
time = range(length)
gamma = 0.99
c0 = 1
# st = np.random.poisson(lam=10, size=length).astype(float)*0.01
st = np.random.normal(loc=0.0, scale=1.0, size=length)*0.05
ct=np.zeros(shape = length)
ctn=np.zeros(shape = length)
spikes = [100, 350, 470]
spikeAmplitude = 1
#自回归模型
for i in range(1, length):
if i in spikes:
ct[i] = gamma*ct[i-1] + spikeAmplitude
else:
ct[i] = gamma*ct[i-1]
ctn[i] = ct[i]+ st[i]
print(ct[i])
fig,ax = plt.subplots()
ax.plot(time, ct, color = 'b', label = 'true')
ax.scatter(time, ctn, color = 'r',s=3, label = 'add noise')
ax.set_xlabel('Frame')
ax.set_ylabel('Calcium Transient')
plt.title('Gamma is {}'.format(gamma))
plt.legend(loc = "best")
fig.show()
--------------------
matalb
len = 500;
gamma = 0.98;% 0.99,200帧;
spikes = [100];
c0 = 1;
ct = zeros(1,len,'double');
ctn = ct;% noise
spikeAmplitude = 2;
mu = 0;sigma = 1;
st = normrnd(mu, sigma, 1, len) .* 0.05;
for i = 2:len
% do something
if ismember(i,spikes)
ct(1,i) = gamma*ct(1,i-1) + spikeAmplitude;
else
ct(1,i) = gamma*ct(1,i-1);
end
ctn(1,i) = ct(1,i)+ st(1,i);
end
figure;
plot(ct);hold on;
plot(ctn);