Fitting exponential decay in Python3

[CharlieNewey]

Hi there,

I have a set of coordinates (data points) that I want to use Python3 to fit an exponential decay curve to. I've used this resource here as a base for building my program.

The problem is, no matter what the x-value I put in is, the y-value ALWAYS comes up as 1.0!

My code is below.

# curve fitting algorithm
# least squares fit method
import math

n = 800
p0 = (5*(10**6))
decay = (1.16*(10**-3))

def calculate_a_fit():
    numerator = 0
    numeratorPos = 0
    numeratorNeg = 0
    denominator = 0
    denominatorPos = 0
    denominatorNeg = 0

    # calculate top half of sum/algorithm
    
    sigma = 0
    for i in range(0, n):        
        y = (p0 * math.e)**(-(decay * n))

        sigma += ((i**2)*y)
    numeratorPos += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (y*(math.log(y)))
    numeratorPos *= sigma
    
    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (i*y)
    numeratorNeg += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (i*y*(math.log(y)))
    numeratorNeg *= sigma

    numerator = numeratorPos - numeratorNeg

    # calculate bottom half of sum/algorithm

    sigma = 0
    for i in range(0, n):        
        y = (p0 * math.e)**(-(decay * n))

        sigma += (y)
    denominatorPos += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += ((i**2)*y)
    denominatorPos *= sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += ((i*y)**2)
    denominatorNeg += sigma

    denominator = denominatorPos - denominatorNeg

    a = (numerator/denominator)

    print(a)
    return a

calculate_a_fit()

def calculate_b_fit():
    numerator = 0
    numeratorPos = 0
    numeratorNeg = 0
    denominator = 0
    denominatorPos = 0
    denominatorNeg = 0

    # calculate top half of sum/algorithm
    
    sigma = 0
    for i in range(0, n):        
        y = (p0 * math.e)**(-(decay * n))

        sigma += (y)
    numeratorPos += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (i*y*(math.log(y)))
    numeratorPos *= sigma
    
    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (i*y)
    numeratorNeg += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += (y*(math.log(y)))
    numeratorNeg *= sigma

    numerator = numeratorPos - numeratorNeg

    # calculate bottom half of sum/algorithm

    sigma = 0
    for i in range(0, n):        
        y = (p0 * math.e)**(-(decay * n))

        sigma += (y)
    denominatorPos += sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += ((i**2)*y)
    denominatorPos *= sigma

    sigma = 0
    for i in range(0, n):
        y = (p0 * math.e)**(-(decay * n))

        sigma += ((i*y)**2)
    denominatorNeg += sigma

    denominator = denominatorPos - denominatorNeg

    b = (numerator/denominator)
    
    print(b)
    return b

a = calculate_a_fit()
b = calculate_b_fit()

# maybe y = (a ** (b * 800))?

for i in range(0, n):
    y = (math.exp(a) * math.e) ** (b * n)
    print (y)

I used a for loop for each sigma.
What am I doing wrong?

Thanks in advance :)

[winmic]

the problem is on line # 157, you are looping from 0 to n, but are using n in the equation. So 'y = (math.exp(a) * math.e) ** (b * n)' should be 'y = (math.exp(a) * math.e) ** (b * i)'(if I understand the math correctly).

you may also want to look at the decimal module if you need to have more exact numbers.

[woooee]

Also, calculate (p0 * math.e) once before the for() instead of every pass calculating it again. You should be using one for() loop, since they all do the same thing.

numeratorPos1 = 0
numeratorPos2 = 0
denominatorPos = 0
p0_math_e = p0 * math.e
for x in range(0, n):        
    y = p0_math_e**(-(decay * x))

    numeratorPos += ((x**2)*y)
    numeratorPos += (y*(math.log(y)))
    denominatorPos += y
#
# etc.

[CharlieNewey]

Hey,

Tried the solution in post #2, although it may have been correct, it didn't properly work; I still have the same problem.

And the idea in post #3 was good, I used the separate for loops to keep it easier for me to understand - a separate for loop for each sigma in the equation - I'll give it a try though, thanks.

Anybody have any ideas?

[LoveMyPadres]

Only other way I know to debug this is to break it down. Check each variable that is part of the final calculation and make sure it makes sense. Mostly likely, there is a subtle bug in how you are doing one of the calculations. That's how I find my bugs anyway. Happy Hunting!

[pyTony]

It is considered bad, because this time for example you had same bug repeated 14 times in your code, for example. Instead insert print functions and any input statement for pause or use debugger break points and watches.

[Gribouillis]

Here is another way to code this

import math

def fsum(f, g, pairs):
    return sum(f(x)*g(y) for (x,y) in pairs)

def one(x):
    return 1.0

def ident(x):
    return x

def square(x):
    return x*x

def plog(x):
    return x * math.log(x)

def fit(pairs):
    pairs = list(pairs)
    si, op, ii, ip, oi = [ fsum(f, g, pairs) for (f, g) in (
            (square, ident),
            (one, plog),
            (ident, ident),
            (ident, plog),
            (one, ident),
        )]
    den = float(oi * si - ii ** 2)
    a = (si * op - ii * ip)/den
    b = (oi * ip - ii * op)/den
    return a, b

def gen_pairs():
    n = 800
    p0 = (5*(10**6))
    decay = (1.16*(10**-3))
    print("expected b: %.16f" % (-decay * math.log(p0*math.e)))
    
    for i in range(n):
        y =  (p0 * math.e)**(-(decay * i))
        yield (i, y)
        
def main():
    a, b = fit(gen_pairs())
    print(a)
    print(b)
    
if __name__ == "__main__":
    main()

"""
my output -->
expected b: -0.0190529402256621
-2.40854223893e-16 # expected a is 0.0
-0.0190529402257
"""

[CharlieNewey]

I found the problem and fixed it - I'd used n in the for loops instead of i. Thanks for all the suggestions on how to clean up my code though, the tips and things are proving useful. Thanks :)