Цей блог представляє рішення вправи 4, глави 1 з книги Теоретична нейронаука Пітера Деяна та Л.Ф. Абботта, з використанням Python.
Я ще не писав жодного блогу для попередніх вправ з цієї книги, тому що вважаю, що ця вправа охоплює всю необхідну теорію та інструменти для вирішення попередніх.
Увесь код, написаний нижче, можна знайти в наступному Google Colab нотатнику. Не соромтесь використовувати його:
[
Google Colab
Опис редагування
colab.research.google.com
](https://colab.research.google.com/drive/14UNbLBV-SGpeckCy18-mzZYOxCo3qrFF?usp=sharing&source=post_page-----8c57b8a2902f--------------------------------)
from sympy.stats import Uniform, sample
from sympy import Symbol, lambdify, log, exp, cos, pi
from sympy.functions.elementary.piecewise import Piecewise, ExprCondPair
import numpy as np
import matplotlib.pyplot as plt
Це всі імпорти, які нам знадобляться для решти блогу.
Отже, без зайвих слів, давайте перейдемо до реалізації!
Аналіз проблеми
Задача полягає ось в чому:
Скріншот з PDF вправ книги
Коротко кажучи, нам потрібно наблизити функцію частоти виведення за допомогою розв'язку диференціального рівняння, але з невеликою зміною: кожного разу, коли відбувається спайк, ми повинні додавати 1/τ до поточної значення частоти виведення. В кінці ми повинні знайти оптимальне значення τ для наближення функції частоти з середньоквадратичною помилкою наближення.
Ми розв’яжемо диференціальне рівняння пізніше. А зараз зосередимось на генерації послідовності спайків.
Генератор спайків Пуассона
Ми можемо використати генератор спайків Пуассона для створення послідовностей спайків з постійною та часозалежною частотою виведення.
Задача полягає в наближенні наступного рівняння:
Де r(t) — це функція частоти виведення, а t — час у мілісекундах. Зверніть увагу, що частота виведення буде коливатися з часом, тому ми будемо робити спостереження за часом залежним від спайків.
Для цього я використаю один з методів, описаних у книзі, хоча є й інші способи досягти того ж результату.
Метод книги розподілений на два основні етапи: генерація послідовності спайків з постійною частотою виведення, що буде визначена як максимальне можливе значення для функції частоти виведення, що залежить від часу. Потім, ця послідовність спайків буде редукована з використанням функції частоти виведення.
_ cos(t)_ коливається між -1 і 1. У часі 0, cos(0) = 1, що є максимальним значенням косинуса.
Це означає, що максимальна частота виведення становить 200 Гц.
З цією інформацією тепер ми можемо запустити наш процес Пуассона.
Процес генерації нових спайків виглядає наступним чином:
Де ti_ — це час спайка i, xrand_ — випадкова змінна з рівномірного безперервного розподілу, а r — максимальне значення функції частоти виведення (200 Гц).
Коротко кажучи, що ми робимо — це генеруємо випадкові інтервали між спайками і додаємо їх до часу останніх спровокованих спайків для генерації нових спайків.
Наступний код генерує і будує графік послідовності спайків з постійною частотою виведення:
#Генерація спайків
a = 0
b = 1
X = Uniform('x', a, b)
spike_train = np.array([0])
duration = 1 #секунда
size_sample = 200
r = 200 #Гц
while True:
samples = -np.log(sample(X, size = size_sample, library = 'numpy'))/r
samples = np.cumsum(np.insert(samples, 0 , spike_train[-1]))
spike_train = np.concatenate((spike_train, samples[1:]))
if spike_train[-1] > duration:
break
spike_train = spike_train[spike_train < duration]
spike_train = spike_train[1:] #Видаляємо перший спайк зі значенням 0
#Побудова графіка спайків
plt.figure(figsize = (20, 5))
plt.plot(spike_train, np.zeros(spike_train.shape[0]), '.k')
Послідовність спайків з постійною частотою виведення
Що робить цей код? Він генерує 200 зразків з рівномірного розподілу на кожному циклі, обчислює log() і ділить на максимальну частоту виведення. Потім починає додавати всі отримані значення для обчислення часу спайків. Після цього перевіряє, чи не перевищено час випробування (1 секунда), і в разі перевищення припиняє виконання циклів. В кінці видаляються всі непотрібні часи спайків.
Щоб перетворити цю послідовність спайків на таку, що залежить від часу, нам потрібно пройти по всіх раніше згенерованих спайках, поділити кожен час спайка, обчислений на основі функції частоти, на максимальну частоту виведення, отримати новий зразок з випадкової змінної і порівняти два значення. Якщо зразок більше за значення з ділення, ми його відкидаємо, інакше — залишаємо.
Цей процес можна виразити наступним рівнянням:
Де ti_ — це час спайка i. Якщо ви заплуталися, я настійно рекомендую звернутися до книги.
Ми можемо реалізувати цей процес за допомогою наступного коду:
def time_dependant_trial(spike):
global r #Максимальне значення частоти виведення
xrand = sample(X, library='numpy')
value = rt(spike) #обчислення функції частоти виведення
return spike if ((value/r) >= xrand) else 0.0 #Перетворюємо на нуль всі неприпустимі спайки
time_dependant_trial_vectorized = np.vectorize(time_dependant_trial) #Векторизація функції для ітерацій по numpy масиву
Функція timedependanttrial фактично виконує весь процес редукції, описаний вище. Зверніть увагу, що відкинуті часи спайків представлені як нулі, які пізніше будуть відрізані.
Оскільки майже неможливо отримати рівно 0 через те, як працюють безперервні ймовірності, це не усуне перший спайк, що був згенерований.
r = 200 #Гц максимальна частота виведення
t = Symbol('t')
rt = 100*(1+cos(2*pi*t/.3)) #функція частоти виведення, що залежить від часу
rt = lambdify(t, rt) #lamdify для швидшого обчислення
time_dependant_sequence = time_dependant_trial_vectorized(spike_train)
time_dependant_sequence = time_dependant_sequence[time_dependant_sequence != 0] #Відрізаємо значення 0
fig, ax = plt.subplots(figsize = (20, 5))
ax.plot(time_dependant_sequence, np.zeros(time_dependant_sequence.shape[0]), '.k')
Послідовність спайків з частотою виведення, що залежить від часу
Аппроксимація функції частоти виведення
Давайте знову подивимося на диференціальне рівняння, яке нам потрібно розв'язати:
Поділивши і помноживши вираз на τ, raprox_ та dt відповідно, ми можемо тепер інтегрувати обидві сторони і отримати наступне:
К — це стала, яку ми пізніше обчислимо.
Потім, застосувавши експоненти, ми отримуємо таке рівняння:
Тепер у нас є можливість вибрати або позитивну, або негативну сторону функції абсолютного значення. Ми виберемо позитивну, оскільки, фізично, негативна частота виведення не існує.
Ми можемо перезаписати рівняння так:
Зверніть увагу, що ми отримали ще одну сталу величину.
Це не зовсім коректно з математичної точки зору, але ми переозначимо сталу k для отримання наступного результату:
Щоб визначити значення k, ми можемо припустити, що в момент часу 0 аппроксимована частота виведення повинна становити 200 Гц, оскільки це значення реальна функція приймає в момент часу 0. Тому k дорівнює 200.
Але ми ще не завершили. Тепер нам потрібно переозначити рівняння щоразу, коли ми отримуємо спайк, і ми зробимо це через трансляцію функції.
Дозвольте пояснити цей процес на прикладі.
Я використаю desmos.com для підтримки пояснення.
Припустимо, що наша функція частоти виведення виглядає такою:
Графік функції частоти виведення
І щоразу, коли ми отримуємо спайк, ми хочемо додавати 1 до значення частоти виведення.
Також припустимо, що ми отримали спайк у момент часу 1.
Частота виведення в момент часу 1 і частота виведення в той самий момент часу, збільшена на 1 одиницю, виглядають наступним чином:
Щоб визначити, в який момент часу ми отримаємо другий результат, потрібно розв’язати таке рівняння:
Або:
Де ts_ — це час, коли було зафіксовано спайк, у нашому випадку це 1.
Якщо розв’язати це рівняння для ts_ = 1, отримаємо наступний результат:
Це означає, що приблизно в момент часу -0.31 ми отримаємо нове значення частоти виведення.
Але ми хочемо, щоб у час 1 частота виведення була такою ж, як і у час -0.31.
Нам потрібно отримати різницю між часом спайка і часом, коли повинна бути досягнута бажана частота виведення, а потім перемістити функцію вправо на отриману величину, як показано нижче:
Примітка: Для того, щоб перемістити функцію вздовж осі x вправо, потрібно відняти одиниці, на які ми будемо переміщувати функцію, від значення x.
Якщо ми підставимо t = 1 у попередню функцію, ми отримаємо бажане значення:
Ми можемо використати наступний код, щоб відтворити попередню процедуру:
k = 200 #Сталу з розв'язку диференціального рівняння
tau_list = [0.01, .02, .1, .05]
r_aprox_list = []
t = Symbol('t')
#Створимо поелементну функцію для кожного tau
for tau in tau_list:
dif_time = 0
r_aprox_aux = exp(-t/tau)
r_aprox_aux_list = []
r_aprox_aux_list.append(ExprCondPair(k*r_aprox_aux, (t >= 0) & (t < time_dependant_sequence[1]))) #Визначаємо частину функції, де ще не було спайка
time_dependant_sequence= np.append(time_dependant_sequence, duration)
for i in range(1,time_dependant_sequence.shape[0]-1):
spike = time_dependant_sequence[i]
spike_next = time_dependant_sequence[i+1]
aux_time = -tau*log(r_aprox_aux + 1/(tau*k)).subs(t, spike) #Отримуємо новий час для трансляції
dif_time = spike - aux_time #Отримуємо різницю між часом спайка і новим часом
r_aprox_aux = exp(-(t-dif_time)/tau) #Переміщаємо графік
expression = ExprCondPair(k*r_aprox_aux, (t >= spike) & (t < spike_next))
r_aprox_aux_list.append(expression)
r_aprox_list.append(lambdify(t, Piecewise(*r_aprox_aux_list))) #Створюємо поелементну функцію
У цьому фрагменті коду ми створюємо поелементну функцію, використовуючи клас Piecewise з бібліотеки SymPy. ExprCondPair — це клас, що дозволяє нам визначати різні частини нашої поелементної функції, які зберігаються в списку і використовуються в кінці як аргументи для створення нового об'єкта Piecewise. Крім того, функція lambdify дозволяє перевести нашу функцію з SymPy в NumPy для швидшого обчислення.
Результати можна візуалізувати за допомогою наступного коду:
xvalues = np.arange(0, 1, 0.001)
f = lambdify(t,100*(1+cos(2*pi*t/.3)))
fig, ax = plt.subplots(figsize = (20, 5))
for i in range(len(tau_list)):
ax.plot(xvalues, r_aprox_list[i](xvalues), label = str(tau_list[i]))
ax.xaxis.set_label_text('Час (мс)')
ax.yaxis.set_label_text('Частота виведення (Гц)')
ax.set_title('Аппроксимована функція частоти виведення')
ax.plot(xvalues, f(xvalues), label = "Оригінальна функція частоти")
ax.legend()
Отримання середньоквадратичної похибки для значень tau
Щоб визначити, яке значення τ найкраще аппроксимує функцію частоти виведення, ми можемо отримати середньоквадратичну похибку для всіх значень τ наступним чином:
xvalues = np.arange(0, 1, 0.001)
squared_error_list = []
for i in range(len(tau_list)):
squared_error = np.mean((f(xvalues) - r_aprox_list[i](xvalues))**2)
squared_error_list.append(squared_error)
for i in range(len(tau_list)):
print(f"Середньоквадратична похибка для tau = {tau_list[i]} становить {squared_error_list[i]}") #tau = 0.02 дає найкращий результат
Найкращу аппроксимацію було отримано при τ = 0.02.
Висновок
Функції, які ми отримали, хоча й не дуже добре аппроксимують функцію частоти виведення, все ж емітують коливальний характер косинусоїдальної функції.
Однак я вважаю, що цінність цього вправи полягає в основному в розв'язуванні задачі. Як я згадував на початку, це завдання дозволяє вам легше вирішувати всі попередні.
Якщо ви дійшли до цього місця, дякую за увагу. Сподіваюся, цей блог був корисний для вас.
Перекладено з: Theoretical Neuroscience, Exercise 4, Chapter 1: Aproximating the firing rate function