Среднеквадратичное значение функции

Question

Теперь полный код / ​​вопросы

Я хотел бы оценить случайные колебания функции v - поэтому я хотел бы рассчитать среднеквадратичное значение ее:

import numpy as np
import matplotlib.pyplot as plt




def HHmodel(I,length, area):

        v = []
        m = []
        h = []
        z = []
        n = []
        squares = []
        vsquare = (-60)*(-60)
        sumsquares = 0
        rms = []
        a= []
        dt = 0.05
        t = np.linspace(0,100,length)


        #constants
        Cm = area#microFarad
        ENa=50 #miliVolt
        EK=-77  #miliVolt
        El=-54 #miliVolt
        g_Na=120*area #mScm-2
        g_K=36*area #mScm-2
        g_l=0.03*area #mScm-2

        def alphaN(v):
            return 0.01*(v+50)/(1-np.exp(-(v+50)/10))

        def betaN(v):
            return 0.125*np.exp(-(v+60)/80)

        def alphaM(v):
            return 0.1*(v+35)/(1-np.exp(-(v+35)/10))

        def betaM(v):
            return 4.0*np.exp(-0.0556*(v+60))

        def alphaH(v):
            return 0.07*np.exp(-0.05*(v+60))

        def betaH(v):
            return 1/(1+np.exp(-(0.1)*(v+30)))

        #Initialize the voltage and the channels :
        v.append(-60)
        rms.append(1)
        m0 = alphaM(v[0])/(alphaM(v[0])+betaM(v[0]))
        n0 = alphaN(v[0])/(alphaN(v[0])+betaN(v[0]))
        h0 = alphaH(v[0])/(alphaH(v[0])+betaH(v[0]))

        #t.append(0)
        m.append(m0)
        n.append(n0)
        h.append(h0)

        #solving ODE using Euler's method:
        for i in range(1,len(t)):
            m.append(m[i-1] + dt*((alphaM(v[i-1])*(1-m[i-1]))-betaM(v[i-1])*m[i-1]))
            n.append(n[i-1] + dt*((alphaN(v[i-1])*(1-n[i-1]))-betaN(v[i-1])*n[i-1]))
            h.append(h[i-1] + dt*((alphaH(v[i-1])*(1-h[i-1]))-betaH(v[i-1])*h[i-1]))
            gNa = g_Na * h[i-1]*(m[i-1])**3
            gK=g_K*n[i-1]**4
            gl=g_l
            INa = gNa*(v[i-1]-ENa)
            IK = gK*(v[i-1]-EK)
            Il=gl*(v[i-1]-El)
            v.append(v[i-1]+(dt)*((1/Cm)*(I[i-1]-(INa+IK+Il))))
            #v.append(v[i-1]+(dt)*((1/Cm)*(I-(INa+IK+Il))))
        meansquare = np.sqrt((np.square(v).sum()))
        return v,area,meansquare




spikeEvents = []  #timing each spike
length = 1000*5  #the time period

fluctuations = []
output = []



for j in range(1, 10):
    barcode = np.zeros(length)
    noisyI = np.random.normal(0,9,length)
    area = 1.0+0.1*j
    res = HHmodel(noisyI,length,area)
    output.append(res[2])




print('Done.')

цель должна состоять в том, чтобы колебания v каким-то образом возрастали с увеличением размера a - я думал, что здесь среднеквадратичная амплитуда является разумной мерой

BR

edit:

 for i in range(1,len(t)):
            m.append(m[i-1] + dt*((alphaM(v[i-1])*(1-m[i-1]))-betaM(v[i-1])*m[i-1]))
            n.append(n[i-1] + dt*((alphaN(v[i-1])*(1-n[i-1]))-betaN(v[i-1])*n[i-1]))
            h.append(h[i-1] + dt*((alphaH(v[i-1])*(1-h[i-1]))-betaH(v[i-1])*h[i-1]))
            gNa = g_Na * h[i-1]*(m[i-1])**3
            gK=g_K*n[i-1]**4
            gl=g_l
            INa = gNa*(v[i-1]-ENa)
            IK = gK*(v[i-1]-EK)
            Il=gl*(v[i-1]-El)
            v.append(v[i-1]+(dt)*((1/Cm)*(I[i-1]-(INa+IK+Il))))
            z.append(v[i-1]-np.mean(v))
            #v.append(v[i-1]+(dt)*((1/Cm)*(I-(INa+IK+Il))))
        mean = sum(np.square(v))/len(v)
        squared_diffs =[(item-mean)**2 for item in v]
        ms_diff = sum(squared_diffs)/len(squared_diffs)
        rms_diff =np.sqrt(ms_diff)
        return v,area,rms_diff

edit2: График для j в диапазоне (1,10) - синий: значение rms, рассчитанное при редактировании 1, желтый 1 / sqrt (j)

edit3:

График для j в диапазоне (1100) - но «размер» колебаний должен увеличиваться, а не уменьшаться и центрироваться где-то

Answer 1

Несколько незначительных замечаний:

• Итак, в основном ваша "функция" v является дискретной оценкой некоторой функции за один шаг, а не истинной функцией, но это не совсем так. здесь актуально.

• Как указано в комментариях выше, вы должны вычислить v для всех временных шагов и агрегировать возведенные в квадрат значения, затем суммировать их за пределами l oop и нормализовать путем деления на len(v).

• Также неясно, почему при итерации i вы вычисляете v[i], но соответствующее вычисленное вами квадратное значение равно v[i-1] в квадрате. Следует использовать тот же индекс на той же итерации l oop, иначе вы, скорее всего, пропустите элемент.

Я бы сказал, что причина в том, что результат бесполезен, заключается в том, что root -средний квадрат на самом деле никогда не используется для выходных данных функции (в данном случае RMS - это просто некое полезное средство, которое придает дополнительный вес выбросам); скорее RMS обычно используется для ошибки или дисперсии выходов этой функции. Среднеквадратичная ошибка или дисперсия говорит вам, насколько далеко в исходных единицах функции среднее значение функции отличается от среднего значения?). Обратите внимание, что это действительно важный показатель c, только если вы ожидаете, что значение v будет постоянным.

Учитывая все это, трудно сказать из вашего вопроса, каково ваше намерение и что вы на самом деле пытаясь сделать с этой информацией, так что я буду догадываться, что вас действительно волнует, насколько значение v отличается от среднего. В этом случае вы можете использовать среднеквадратичное отклонение от среднего значения v, рассчитанного так:

for i in range(1,len(t)):
        #calculate v[i] here, omitted for simplicity

    # get mean value
    mean = sum(squares)/len(squares)

    # you want to get the squared value of the difference, not the value itself
    squared_diffs = [(item - mean)**2 for item in v)]


    # get mean squared diff
    ms_diff = sum(squared_diffs) / len(squared_diffs)

    # return root of mean squared diff
    rms_diff = np.sqrt(ms_diff)

    return v,area,rms_diff

Опять же, это полезно, только если вы ожидаете, что выходные данные v будут постоянными. Если нет, вы попытаетесь подогнать к функции другую модель (линейную, квадратичную c и др. c.), А затем рассчитать среднеквадратичную ошибку. Вопрос был бы гораздо понятнее, если бы вы указали цель этого расчета.