0

I asked a similar question in January that @Miłosz Wieczór was kind enough to answer. Now, I am faced with a similar but different challenge since I need to fit two parameters (fc and alpha) simultaneously on two datasets (e_exp and iq_exp). I basically need to find the values of fc and alpha that are the best fits to both data e_exp and iq_exp.

import numpy as np
import math
from scipy.optimize import curve_fit, least_squares, minimize

f_exp  = np.array([1, 1.6, 2.7, 4.4, 7.3, 12, 20, 32, 56, 88, 144, 250000])
e_exp  = np.array([7.15, 7.30, 7.20, 7.25, 7.26, 7.28, 7.32, 7.25, 7.35, 7.34, 7.37, 11.55])
iq_exp = np.array([0.010, 0.009, 0.011, 0.011, 0.010, 0.012, 0.019, 0.027, 0.038, 0.044, 0.052, 0.005])

ezero  = np.min(e_exp)
einf   = np.max(e_exp)

ig_fc     = 500
ig_alpha  = 0.35

def CCRI(f_exp, fc, alpha):
    x   = np.log(f_exp/fc)
    R  = ezero + 1/2 * (einf - ezero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (einf - ezero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    RI = np.sqrt(R ** 2 + I ** 2)
    return RI

def CCiQ(f_exp, fc, alpha):
    x   = np.log(f_exp/fc)
    R  = ezero + 1/2 * (einf - ezero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (einf - ezero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    iQ = I / R
    return iQ

poptRI, pcovRI = curve_fit(CCRI, f_exp, e_exp, p0=(ig_fc, ig_alpha))

poptiQ, pcoviQ = curve_fit(CCiQ, f_exp, iq_exp, p0=(ig_fc, ig_alpha))

einf, ezero, and f_exp are all constant plus the variables I need to optimize are ig_fc and ig_alpha, in which ig stands for initial guess. In the code above, I get two different fc and alpha values because I solve them independently. I need however to solve them simultaneously so that fc and alpha are universal.

  • Is there a way to solve two different functions to provide universal solutions for fc and alpha?
naughty_waves
  • 265
  • 1
  • 13
  • http://scipy-lectures.org/intro/scipy/auto_examples/plot_2d_minimization.html create a objecftive function that uses both results from your curve_fits ?? – Joe Apr 15 '20 at 11:53
  • Thank you for commenting and providing a link, @Joe. I am still not sure how I go about to actually do it in terms of combining the two into an objective function. Could you perhaps elaborate? – naughty_waves Apr 15 '20 at 16:31

1 Answers1

1

The docs state on the second returned value from curve_fit:

pcov

The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov)).

So if you want to minimize the overall error, you need to combine the errors of both your fits.

def objective(what, ever):

    poptRI, pcovRI = curve_fit(CCRI, f_exp, e_exp, p0=(ig_fc, ig_alpha))

    poptiQ, pcoviQ = curve_fit(CCiQ, f_exp, iq_exp, p0=(ig_fc, ig_alpha))

    # not sure if this the correct equation, but you can start with it
    err_total = np.sum(np.sqrt(np.diag(pcovRI))) + np.sum(np.sqrt(np.diag(pcoviQ)))

    return err_total

On total errors of 2d Gaussian functions:

https://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/

Update: Since you want poptRI and poptiQ to be the same, you need to minimize their distance.

This can be done like

from numpy import linalg

def objective(what, ever):

    poptRI, pcovRI = curve_fit(CCRI, f_exp, e_exp, p0=(ig_fc, ig_alpha))

    poptiQ, pcoviQ = curve_fit(CCiQ, f_exp, iq_exp, p0=(ig_fc, ig_alpha))

    delta = linalg.norm(poptiQ - poptRI)
    return delta

Minimizing this function will (should) result in similar values for poptRI and poptiQ. You take the parameters as vectors, and try to minimize the length of their delta vector.

However, this approach assumes that poptRI and poptiQ (and their coefficients) are about in the same range since you are using some metric on them. If say one if them is in the range 2000 and the other in the range 2. Then the optimizer will favour tuning the first one. But maybe this is fine.

If you somehow want to treat them the same you need to normalize them. One approach (assuming all coefficients are similar) could be

linalg.norm((poptiQ / linalg.norm(poptiQ)) - (poptRI / linalg.norm(poptRI))))

You normalize the results to unit vectors, then subtract them, then create the norm.

The same is true for the inputs to the function, but it might not be that important there. See the links below.

But this strongly depends on the problem you are trying to solve. There is no general solution.

Some links related to this:

Is normalization useful/necessary in optimization?

Why do we have to normalize the input for an artificial neural network?

Another objective function:

It this what you are trying to do? You want to find the best fc and alpha so the fit results of both functions are as close as possible?

def objective(fc, alpha):

    poptRI, pcovRI = curve_fit(CCRI, f_exp, e_exp, p0=(fc, alpha))

    poptiQ, pcoviQ = curve_fit(CCiQ, f_exp, iq_exp, p0=(fc, alpha))

    delta = linalg.norm(poptiQ - poptRI)
    return delta
Joe
  • 6,758
  • 2
  • 26
  • 47
  • Thank you again, @Joe. I actually want to find the best estimates for `fc` and `alpha` which I guess occur when the errors are the lowest possible values. So, instead of combining `pcovRI` and `pcoviQ` I need to find the common value comprising both `poptRI` and `poptiQ`. Do you have an idea on how to combine them? – naughty_waves Apr 15 '20 at 19:40
  • See the second objection function. – Joe Apr 15 '20 at 20:54
  • Thank you, @Joe. You are probably two steps (if not more) ahead of me. If run them independently `poptRI[0]=2549.59`, `poptRI[1]=0.2921`, `poptiQ[0]=2801.64`, and `poptiQ[1]=0.3202` in which `[0]=fc` and `[1]=alpha`. If I run your objective function with `what=e_exp` and `ever=iq_exp`, delta returns `252.04` which I expected it to be somewhere between `poptRI[0]` and `poptiQ[0]`. I also need `alpha` (`poptRI[1]` and `poptiQ[1]`) but this is not provided. I tried to simply add and return `delta2=linalg.norm(poptiQ[1] - poptRI[1])` but this value is also too small. I apologize for my ignorance. – naughty_waves Apr 16 '20 at 09:48
  • I did not completely understand what you were trying to do. You might have to add more free variables to the objective function like `def objective(what, ever, you, need)`. Maybe `alpha` and `fc`? I also added something on normalization to the question. – Joe Apr 16 '20 at 10:10
  • So, since `fc`and `alpha` from `CCRI` (`poptRI[0] ` and `poptRI[1]`) are different to the same parameters from `CCiQ` (`poptiQ[0] ` and `poptiQ[1]`), I would like to optimize both of them simultaneously in order to find the best universal fits for `fc`and `alpha` that I can put into my model. Is `delta` the normalized `alpha`? If that is the case, how do I back-calculate to find the real value of `alpha`? What about `fc`? Again, I am sorry for being annoyingly ignorant. I am just trying to wrap my thick head around it. I feel bad for wasting your time because I am clearly not worthy of it. – naughty_waves Apr 16 '20 at 17:43
  • It's completely fine, you are not wasting my time. Basically you just need to define what you would like to optimize and which parameters define the result. You could for example use `def objective(fc1, fc2, alpha1, alpha2)`. As mentioned above, only you fully understand what you want to do. Another option would be to use an `inner` and an `outer` optimization, basically a `minimize` which another `minimize` in its objective function. – Joe Apr 16 '20 at 18:16
  • About how to get the parameters, you called it `back-calculate`, you could e.g. use a callback function (see keyword `callback` in `minimize`) to be called after each iteration and write the results (parameters + value of objective function) to a file. Then you can track how your optimization converges (hopefully) and you have the parameters. Or if you are looking for some inner variable used during the call of the objective function write the values to a file. Or append to a list or an array. – Joe Apr 16 '20 at 18:22
  • Thank you, @Joe. `delta` is basically the distance between `poptRI[0]` and `poptIQ[0]`. Instead of the distance between the two, I want to find universal `popt[0]` and `popt[1]` (`fc` and `alpha`) that applies to both since I can only use singular `fc` and `alpha` in my model. In my thick head, there should be a `popt[0]` and `popt[1]` that is somewhat in between the ones I get when I `curve_fit` them independently (like I did in the beginning). Is it even possible to fit a curve that depends on two variables (`fc` and `alpha`) to two different curves? I realize I am terrible at explaining. – naughty_waves Apr 17 '20 at 16:13
  • Given these two functions, use, for example, `curve_fit` to find the optimal `fc` and `alpha` that applies to both curves. Then, use `fc` and `alpha` to fit the theoretical model to your experimental data. Ideally, I would throw both `CCRI` and `CCiQ` into `curve_fit` which would directly provide the optimized `fc` and `alpha` but that is alas not possible. – naughty_waves Apr 17 '20 at 16:21
  • I added another objective function at the bottom of my answer. That what you are looking for? – Joe Apr 17 '20 at 16:55
  • Another thing that you could do if you need to somehow dynamically adjust the functions is to define them in the objective function and to define variables needed by them also in there so they are the same for both functions. – Joe Apr 17 '20 at 17:31
  • Thank you for still showing interest, @Joe. The last `objective` function provides the same answer as the other `objective` functions. They provide the distance between `poptRI[0]` and `poptiQ[0]` but not the best fits for `fc` and `alpha` that I need to insert in my model. I would like to define a function that fits the best possible `fc` and `alpha` for both for `e_exp` and `iq_exp`. For example, the averages of `poptRI[0]` and `poptiQ[0]` as well as `poptRI[1]` and `poptIQ[1]` provides a decent fit. This is not really scientifically, and there must be a better way to do it. – naughty_waves Apr 21 '20 at 15:16
  • Perhaps this is what you meant about dynamically adjusting the functions. Again, I apologize for not being on the same page as you. As I said, I think you are a couple of steps ahead of me. I am not very good at this. – naughty_waves Apr 21 '20 at 15:19
  • The last objective function optimizes `fc` and `alpha`. You need to define what makes it a good solution. e.g. the last objective function needs the results of the fit to be similar. Another option as in the first objective function is that both fits are as good as possible, means: work well for both functions. – Joe Apr 21 '20 at 16:27
  • But you can think of other options. e.g. you could use the overall minimum of the parameters for both functions. Or you could flip the sign to maximize something, e.g. say you want to maximize `poptiQ + poptRI`, then since there are only minimization algorithms you would make the objective function return `-(poptiQ + poptRI)`, or wrap that each of them additionally in `np.abs`. What ever you like. So, what defines the ideal `fc` and `alpha` combination for you? – Joe Apr 21 '20 at 16:32