Fehler zu Realdaten minimieren

[claus_kurt]

Liebe Leute,

ich habe eine Funktion

Code: Alles auswählen

def model(x,y):
    #some calculations
    return count

count gibt die Anzahl der aufgetretenen Ereignisse wieder, die in model(x,y) durch 4 Differentialgleichungen bestimmt werden.

Nun kenne ich durch ein anderes Modell simulierte Outputs also ich weiß zb, dass

model(1,0.92*y) = 28
model(1,0.95*y) = 37
model(1,1*y) = 50
...

nun möchte ich gerne in meinem Modell den Parameter y so finden, dass der Fehler zu den Daten klein wird - also y optimal finden.

Dazu habe ich folgendes gemacht

Code: Alles auswählen

root = 38.69140625
xdata = np.array([0.92*root,root,1.05*root])
ydata = np.array([28,100,110])
y0=np.array([1,1,1])

def g(y,xdata,ydata):
    return ydata - HH_model(xdata,y)


fit = optimize.leastsq(g, y0, args=(xdata, ydata))

aber ich bekomme als error

ValueError: could not broadcast input array from shape (3) into shape (1)

Was ist das problem?

Danke und LG

[Sirius3]

Weiter lautete die komplette Fehlermeldung und was ist HH_model?

[claus_kurt]

Code: Alles auswählen

Traceback (most recent call last):

  File "...", line 172, in <module>
    fit = optimize.leastsq(g, y0, args=(xdata, ydata))

   File "...", line 388, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)

   File "..."", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))

   File "..."", line 169, in g
    return ydata - HH_model(xdata,y)

  File "..."", line 125, in HH_model
    v[i + 1] = vT + ((-IIon + IStim) / C) * dt   # in ((uA / cm ^ 2) / (uF / cm ^ 2)) * ms == mV

ValueError: could not broadcast input array from shape (3) into shape (1)

HH_model bestimmt in Abhängigkeit zweier Inputdaten die Spannung v und retourniert den Wert wie oft diese Spannung v über einem bestimmten Grenzwert liegt.

[Sirius3]

Und genau in dieser Funktion tritt auch der Fehler auf, also mußt Du auch dort suchen, ob die Dimensionen der Input-Variablen Deinen Erwartungen entsprechen.
Dass da ein i+1 vorkommt, sieht komisch aus.

[claus_kurt]

v hängt natürlich von der Zeit hab - also die Spannung wird in Abhängigkeit der Zeit berechnet und dann geschaut, ob zu irgendeinem Zeitpunkt v(t) > grenze ist.
Wenn ja, dann wird count um 1 erhöht und beispielsweise bei 50 runs, kommt dann count = 30 raus und das ist der return wert dieser Funktion.

[Sirius3]

Nochmal deutlicher: die Modellfunktion hat einen Fehler, die Modellfunktion kennen wir hier aber nicht, ergo können wir nicht helfen.

[claus_kurt]

Code: Alles auswählen

def HH_model(I,area_factor):
    
    count = 0 
    
    for j in range(0,runs):
        
        # Temporal and stimulus parameters
        t_end = 10  # in ms
        tDel = 1     # in ms
        tDur = 0.3    # in ms
        dt = 0.01   # in ms
        I = I 
        area_factor = area_factor
        
      
        kT = 12             # Temp. factor
        
        #geometry
        d = 2       # diameter in um
        r = d/2     # Radius in um
        l = 10      # Length of the compartment in  um
        A = (2 * np.pi * r * l * 1e-8)*area_factor          # surface [cm^2]
        C = c * A    # uF
        ### ------------- Step 2: COMPUTE ADDITIONAL PARAMETERS -------------
        # Introduction functions of the channel dynamics and Initial valus
        
        
        def alphaM(v): return kT * ((2.5 - 0.1 * (v)) / (np.exp(2.5 - 0.1 * (v)) - 1))
        
        
        def betaM(v):  return kT * (4 * np.exp(-(v) / 18))
        
        
        
        def betaH(v): return kT * (1 / (np.exp(3 - 0.1 * (v)) + 1))
        
        
        def alphaH(v): return kT * (0.07 * np.exp(-(v) / 20))
        
        
        def alphaN(v): return kT * ((1 - 0.1 * (v)) / (10 * (np.exp(1 - 0.1 * (v)) - 1)))
        
        
        def betaN(v): return kT * (0.125 * np.exp(-(v) / 80))
        
        
        # Compute time step
        tSteps = t_end/dt+1
        timeStep = np.arange(0, t_end + dt, dt)
        
        # Compute initial values
        v0 = 0
        m0 = alphaM(v0)/(alphaM(v0)+betaM(v0))
        h0 = alphaH(v0)/(alphaH(v0)+betaH(v0))
        n0 = alphaN(v0)/(alphaN(v0)+betaN(v0))
        
        # Allocate memory for v, m, h, n
        v = np.zeros((int(tSteps), 1))
        m = np.zeros((int(tSteps), 1))
        h = np.zeros((int(tSteps), 1))
        n = np.zeros((int(tSteps), 1))
        
        # Set Initial values
        v[:, 0] = v0
        m[:, 0] = m0
        h[:, 0] = h0
        n[:, 0] = n0
         
        
        ### Noise component
        knoise=  0.003  #uA/(mS)^1/2
        ###  --------- Step3: SOLVE
        for i in range(0, int(tSteps)-1, 1):
        
        # Get current states
           vT = v[i]
           mT = m[i]
           hT = h[i]
           nT = n[i]
        
        # Stimulus current
           IStim = 0
           if tDel / dt <= i <= (tDel + tDur) / dt:
               IStim = I * A    # in uA
           else:
               IStim = 0
        
        
        #  Compute change of m, h and n 
               m[i + 1] = (mT + dt * alphaM(vT)) / (1 + dt * (alphaM(vT) + betaM(vT)))
               h[i + 1] = (hT + dt * alphaH(vT)) / (1 + dt * (alphaH(vT) + betaH(vT)))
               n[i + 1] = (nT + dt * alphaN(vT)) / (1 + dt * (alphaN(vT) + betaN(vT)))
        
        
        # Ionic currents
           iNa = gNa * m[i + 1] ** 3. * h[i + 1] * (vT - vNa)    # in (mS / cm ^ 2) * mV == uA / cm ^ 2
           iK = gK * n[i + 1] ** 4. * (vT - vK)                     # in (mS / cm ^ 2) * mV == uA / cm ^ 2
           iL = gL * (vT-vL)                                           # in (mS / cm ^ 2) * mV == uA / cm ^ 2
           Inoise = (np.random.normal(0, 1) * knoise * np.sqrt(gNa * A))  # in uA
           IIon = ((iNa + iK + iL) * A) + Inoise   # in uA
        # Compute change of v
           v[i + 1] = vT + ((-IIon + IStim) / C) * dt   # in ((uA / cm ^ 2) / (uF / cm ^ 2)) * ms == mV
        
        
        # adjust to resting voltage
        v = v + vRest
     
        # test if there was a spike
        
        if max(v[:]-vRest) > 60:
            count += 1
        
              
           
    return count

[claus_kurt]

Ich habe die Funktion nun angehängt.

[Sirius3]

Konstanten und Unterfunktionen definiert man nicht innerhalb einer Funktion innerhalb einer Schleife. Einige Konstanten haben Namen bekommen, andere nicht, warum?
Einige Einrückungen sind auch nicht korrekt. Die Kommentare sind auch falsch eingerückt, so dass man sich ständig fragen muß, was ist nun korrekt.
So eine lange Funktion schreibt man doch nicht an einem Stück herunter und hofft, dass sie das richtige Ergebnis liefert. Sowas entwickelt man Stück für Stück und testet jeden Schritt ausgiebig.
I und area_factor sind Vektoren. Damit auch A und C, und IStim und Inouse und IIon, v[i+1] aber nicht.

[claus_kurt]

Die Funktion liefert das richtige Ergebnis, aber I und area_factor sind eigentlich nur Zahlen - für jede Zahl wird die Funktion entsprechend der Anzahl von runs durchlaufen.
runs = 30 :
Also I = 10 und area factor = 1
dann wird die Funktion für dieses I und diesen area_factor 30mal ausgeführt und am Ende zurückgegeben wie oft die Grenze also max(..) > 60 überschritten wurde.

Nun möchte ich I schrittweise reduzieren I*0.9, I*0.8 etc und kenne dafür die Treffer aus einem anderen Modell. Ich möchte den area_factor nun so wählen, dass die Abweichungen möglichst gering sind.

[Sirius3]

Du DENKST, dass I und area_factor nur Zahlen sind, aber der Computer weiß es besser.

[claus_kurt]

Wie kann ich denn das beheben Sirius? Wie kann ich die "Realdaten" so übergeben, dass in HH_model kein Dimensionsproblem in den Differentialgleichungen auftritt?

[Sirius3]

Was Du genau machen willst, weiß ich ja nicht. v braucht die richtige Dimension und damit wird count auch zu einem Vektor.

[claus_kurt]

Also das Ziel ist area_factor so zu wählen, dass count mit den beobachteten Häufigkeiten möglichst genau übereinstimmt.
ich übergebe also beispielsweise 5 Beobachtungsergebnisse - dabei kenne ich jeweils I, Aber area_factor ist unbekannt. Dieser soll so gewählt werden, dass es einen area_factor gibt, dass für dieses eine I, die Abweichungen von count (also die Anzahl an Schwellwertüberschreitungen) und die beobachtete Anzahl an Überschreitungen (die ich kenne) minimal wird.