0

I want to fit some data to a curve, using this as a cost function:

def cost_func(x):
    return ((unknown_conc-x[1]*(x[0]*conc_A+
           (1-x[0])*conc_B))**2).sum()

It works when using scipy.optimize, but I want to use statsmodels instead. However I'm struggling with defining a statsmodels formula. Do you have any ideas how to do this?

I tried something like this, but it does not work with this x*A + (1-x)*B:

result = sm.ols(formula="A ~ I(B + C) -1", data=df).fit()
phenix
  • 3
  • 2

1 Answers1

1

Statsmodels/patsy formulas are a language for writing linear models, so you need to find a way to phrase your problem as a formula where the predicted value is linear function of the parameters you want to fit.

In this case, you're doing least-squares fitting where the prediction is (using Python syntax):

x[1]*(x[0]*conc_A + (1 - x[0])*conc_B)

Expanding the terms, we get:

x[1]*x[0]*conc_A + x[1]*(1 - x[0])*conc_B

Let's define new parameters param0 = x[1]*x[0] and param1 = x[1]*(1 - x[0]). Now our prediction becomes

param0*conc_A + param1*conc_B

Notice that these are invertible, i.e. these equalities hold:

x[0] = param0 / (param0 + param1)
x[1] = param0 + param1

so this reparametrization isn't changing the underlying model we're fitting, it's just changing how we represent it. But the new representation is linear in the parameters, so now we can convert it to a statsmodels/patsy formula:

"conc_A + conc_B - 1"

And finally, let's put the value we're fitting our predictions against in the formula, giving:

result = sm.ols("unknown_conc ~ conc_A + conc_B - 1", data=df).fit()

If you fit this you'll get values for param0 and param1, and if you use the equations above you can convert these back to x[0] and x[1] values to compare to what you were getting before.

Nathaniel J. Smith
  • 11,613
  • 4
  • 41
  • 49
  • Thank you, that works just fine and gives me the same results as with scipy. However, do you know how to deal with the variances of the params? So how can I calculate the variances for x[0] and x[1] from the variances for param0 and param1? – phenix Feb 19 '18 at 13:37
  • That's a different and much deeper question... I believe that if you make some of the standard distributional assumptions then OLS statistics will give you an estimate for x[1]'s variance, but not x[0]. Maybe bootstrapping? You should probably ask that on the stats stackexchange. – Nathaniel J. Smith Feb 19 '18 at 18:25