函数的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振幅是一种合理的测量方法

溴

编辑：

 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中计算的rms值，黄色1/sqrt（j）

编辑3:

在范围（1100）内绘制j图-但是波动的“大小”应该增加，而不是减少并集中在某个地方

Answer 1

一些小提示：

• 所以，基本上你的“函数”v是对某个函数的一次离散求值，而不是一个真正的函数，但这在这里并不重要

• 如上面的注释所示，您应该计算所有时间步的v，并聚合平方值，然后在循环之外对它们求和，并通过除以len(v)进行规范化

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

我想说的是，结果没有用处的原因是，均方根实际上从未用于函数的输出（在这种情况下，RMS只是某种不太有用的平均值，为异常值提供了额外的权重）；相反，RMS通常用于该函数输出的误差或方差。RMS error（均方根误差）或variance（方差）告诉您，在函数的原始单位中，平均函数值与平均值的差异有多远？）。请注意，如果您希望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误差。若你们指出了这个计算的目标，那个么这个问题会更清楚