3

I want to estimate the parameters of an RL model from a behavioral dataset with pyomo.

#dummy data
dis_data =  pd.DataFrame([0,1,0,0,0,1], columns=['reward'])
dis_data['Expt']=str(1)
dis_data = dis_data.set_index('Expt')


def max_likelihood_four(x,data):
    lr=x.lr
    sigma=x.sigma
    time = data.shape[0]
    v = np.zeros(([2, time]))
    pr = np.zeros(([2, time]))+0.5
    pr_log = np.zeros(([time]))

    for t in range(time-1):

        pr[1, t] =  1 / (1 + np.exp(-(v[1, t] - v[0, t])/ sigma))
        pr[0, t] = 1-pr[1, t]

        outcome=int(data.ix[t,'reward'])
        v[choice, t + 1] = v[choice, t] + lr * (outcome - v[choice, t])
        v[1 - choice, t + 1] = v[1 - choice, t]

        pr_log[t] = np.log(pr[choice, t])
        # print(i)
    return -np.sum(pr_log)
def pyomo_fit1(df):

    mdl = ConcreteModel()
    mdl.lr    = Var(initialize=np.random.random(1), bounds=(0, 1))
    mdl.sigma = Var(initialize=np.random.random(1), bounds=(0, 10))

    
    residuals = max_likelihood_four(mdl,df)

    mdl.obj = Objective(expr=residulas, sense=minimize)
    SolverFactory('ipopt').solve(mdl)
    return [mdl.lr,mdl.sigma]

parameter_values_1, r1 = pyomo_fit1(dis_data)

I get this error:

"Cannot treat the scalar component 'sigma' as an indexed component"

faraa
  • 575
  • 3
  • 14
  • 42
  • 1
    Can you attach the full error message? (With traceback) – Tamir May 01 '21 at 11:49
  • can only guess because not familiar with pyomo and related applications - quick searches for similar 'scalar' vs. 'indexed component' issues seem to relate to numpy error regarding float issues - see following (web-archived) pages: https://web.archive.org/web/20210501150945/https://groups.google.com/g/pyomo-forum/c/ajySJXA1ol4?pli=1 and https://web.archive.org/web/20210501151044/https://github.com/Pyomo/pyomo/issues/31 – JB-007 May 01 '21 at 15:10

1 Answers1

0

np.random.random(1) returns a numpy.ndarray not an integer or float. When you try to initialize mdl.sigma as a scalar pyomo variable, you need to use a scalar value, or in this case np.random.random(1)[0] (the only element in the array).

If you replace your code for mdl.sigma (and also for mdl.lr) with:

mdl.sigma = Var(initialize=np.random.random(1)[0], bounds=(0, 10))

then the error ("Cannot treat the scalar component 'sigma' as an indexed component") error should disappear.

mdl = ConcreteModel()
mdl.lr    = Var(initialize=np.random.random(1)[0], bounds=(0, 1))
mdl.sigma = Var(initialize=np.random.random(1)[0], bounds=(0, 10))
Tom Krumpolc
  • 56
  • 5