函数的 RMS 值

Question

现在是完整的代码/问题

我想估计函数 v 的随机波动 - 因此我想计算它的 RMS 值：

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 的大小以某种方式增加 - 我在这里认为 rms 幅度是一个合理的度量

BR

编辑：

 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

编辑2： 绘制范围 (1,10) 中的 j - 蓝色：在编辑 1 中计算的 rmsvalue，黄色 1/sqrt(j)

edit3：

在范围（1,100）内绘制 j - 但波动的“大小”应该增加，而不是减少并在某处居中

Answer 1

一些小笔记：

• 所以，基本上你的“函数”v 是对某个函数的一次性离散评估，而不是真正的函数，但这在这里并不真正相关。

1234563

• 还不清楚为什么在迭代 i 中您计算 v[i] 但您计算的相应平方值是 v[i-1] 平方。应该在相同的循环迭代中使用相同的索引，否则您可能最终会丢失一个元素。

我会说结果没有用的原因是均方根并没有真正用于函数的输出（在这种情况下，RMS 只是某种不太有用的平均值，它为异常值提供了额外的权重） ;而 RMS 通常用于该函数输出的误差或方差。 RMS 误差或方差告诉您在函数的原始单位中，平均函数值与平均值相差多远？）。请注意，如果您希望 v 的值保持不变，这才是真正重要的指标。

鉴于这一切，很难从你的问题中说出你的意图是什么以及你实际上试图用这些信息做什么，所以我猜你真正关心的是 v 的值与均值。在这种情况下，您可以使用与 v 平均值的 RMS 差值，计算如下：

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 的输出为常数时才有用。如果不是，您将尝试为函数拟合不同的模型（线性、二次等），然后计算 RMS 误差。如果您指出此计算的目标，问题会更清楚。