3

I am attempting to extract Weibull distribution parameters (shape 'k' and scale 'lambda') that satisfy a certain mean and variance. In this example, the mean is 4 and the variance is 8. It is a 2-unknowns and 2-equations type of problem.

Since this algorithm works with Excel 2010's GRG Solver, I am certain it is about the way I am framing the problem, or potentially, the libraries I am using. I am not overly familiar with optimization libraries, so please let me know where the error is.

Below is the script:

from scipy.optimize import fmin_cg
import math

def weibull_mu(k, lmda):                  #Formula can be found on wikipedia
    return lmda*math.gamma(1+1/k)
def weibull_var(k, lmda):                 #Formula can be found on wikipedia
    return lmda**2*math.gamma(1+2/k)-weibull_mu(k, lmda)**2

def min_function(arggs):
    actual_mean = 4                          # specific to this example
    actual_var = 8                           # specific to this example
    k = arggs[0]
    lmda = arggs[1]
    output = [weibull_mu(k, lmda)-(var_wei)]
    output.append(weibull_var(k, lmda)-(actual_var)**2-(actual_mean)**2)
    return output

print fmin(min_function, [1,1])

This script gives me the following error:

[...]
  File "C:\Program Files\Python27\lib\site-packages\scipy\optimize\optimize.py", line 278, in fmin
    fsim[0] = func(x0)
ValueError: setting an array element with a sequence.
Andre Silva
  • 4,782
  • 9
  • 52
  • 65
TimY
  • 5,256
  • 5
  • 44
  • 57

2 Answers2

4

As far as I can tell, min_function returns a multi-dimensional list, but fmin and fmin_cg does expect that the objective function returns a scalar, if I am not mistaken.

If you are searching the root of the two-equations problem, I suppose it is better that you apply the root function instead. As far as I have been able to find out, scipy does not provide any general optimizers for vector functions.

Anders Gustafsson
  • 15,837
  • 8
  • 56
  • 114
2

I managed to get it to work thanks to Anders Gustafsson's comment (thank you). This script now works if one returns only a scalar (in this case I used something along the lines of least-squares). Also, bounds were added by changing the optimization function to "fmin_l_bfgs_b" (again, thanks to Anders Gustafsson).

I only changed the min_function definition relative to the question.

from scipy.optimize import fmin_l_bfgs_b
import math

def weibull_mu(k, lmda):
    return lmda*math.gamma(1+1/k)
def weibull_var(k, lmda):
    return lmda**2*math.gamma(1+2/k)-weibull_mu(k, lmda)**2

def min_function(arggs):
    actual_mean = 4.                    # specific to this example
    actual_var = 8.                     # specific to this example
    k = arggs[0]
    lmda = arggs[1]
    extracted_var = weibull_var(k, lmda)
    extracted_mean = weibull_mu(k, lmda)
    output = (extracted_var - actual_var)**2 + (extracted_mean - actual_mean)**2
    return output

print fmin_l_bfgs_b(min_function, best_guess, approx_grad = True, bounds = [(.0000001,None),(.0000001,None)], disp = False)

Note: Please feel free to use this script for your own or professional use.

TimY
  • 5,256
  • 5
  • 44
  • 57
  • If you want to apply bounds, use for example [fmin_l_bfgs_b](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_l_bfgs_b.html), that allows you to set bounds on your variables. Also, if you found my answer useful, do not hesitate to upvote it :-) – Anders Gustafsson Jun 28 '12 at 17:58
  • I would love to but I apparently need to have 15 reputation. Also, fmin_l_bfgs gives me an error as well: File"C:\Program Files\Python27\lib\site-packages\scipy\optimize\lbfgsb.py", line 150, in func_and_grad f, g = func(x, *args) TypeError: 'numpy.float64' object is not iterable – TimY Jun 29 '12 at 08:15
  • Typical :-) Looking at the argument list, it seems like you need to set the `approx_grad` argument to `True`. If you don't, either `fprime` needs to defined or the gradient has to be included in the return value from `func`. Does this help? – Anders Gustafsson Jun 29 '12 at 08:37
  • That works perfectly! Thank you so much, Andres. (Also, I just upped your answer since my rep went up) – TimY Jun 29 '12 at 09:37