ValueError: too many values to unpack (expected 3)?

Question

I have been having issues with the code I am trying to right with the model I am trying to code the following error has appeared and being a relative novice I am unsure of how to resolve it.

ValueError                                Traceback (most recent call last)
<ipython-input-2-5f21a0ce8185> in <module>()
     26         proposed[j] = proposed[j] + np.random.normal(0,propsigma[j])
     27         if (proposed[j]>0): # automatically reject moves if proposed parameter <=0
---> 28             alpha = np.exp(logistic_loglik(proposed,time,ExRatio,sig)-logistic_loglik(par_out[i-1,],time,ExRatio,sig))
     29             u = np.random.rand()
     30             if (u < alpha):

<ipython-input-2-5f21a0ce8185> in logistic_loglik(params, t, data, sig)
      3 # set up a function to return the log likelihood
      4 def logistic_loglik(params,t,data,sig):
----> 5     return sum(norm.logpdf(logistic(data, t, params),sig))
      6 
      7 # set standard deviations to be 10% of the population values

<ipython-input-1-c9480e66b7ef> in logistic(x, t, params)
      6 
      7 def logistic(x,t,params):
----> 8     S, R, A = x
      9     r, Nmax, delta_s, beta, gamma, delta_r, delta_a, Emax, H, MICs, MICr = params
     10     N = S + R

ValueError: too many values to unpack (expected 3) 

The model I am trying to code is an MCMC to fit some ODE's to some data I have added the code below for context.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

%matplotlib inline

def logistic(x,t,params):    
    S, R, A = x
    r, Nmax, delta_s, beta, gamma, delta_r, delta_a, Emax, H, MICs, MICr = params
    N = S + R
    E_s = 1 - (Emax * A**H)/(MICs**H + A**H)
    E_r = 1- (Emax * A**H)/(MICr**H + A**H)

    derivs = [r * (1 - N / Nmax ) * E_s * S - delta_s * S - ((beta * S * R)/N), 
     r * (1 - gamma) * (1 - N/Nmax) * E_r * R  - delta_r * R + ((beta * S * R)/N), - delta_a * A]


     return derivs

r = 0.5
Nmax = 10**7
delta_s = 0.025
beta = 10**-2
gamma = 0.5
delta_r = 0.025
delta_a = 0.003
Emax = 2
H = 2
MICs = 8
MICr = 2000

[r, Nmax, delta_s, beta, gamma, delta_r, delta_a, Emax, H, MICs, MICr] = params
S = 9 * 10**6
R = 10**5
A = 5.6
x0 = [S, R, A]

maxt = 2000
tstep = 1
t = np.arange(0,maxt,tstep)

def logistic_resid(params,t,data):
    return logistic(params,t)-data

logistic_out = odeint(logistic, x0, t, args=(params,))

time = np.array([0, 168, 336, 504, 672, 840, 1008, 1176, 1344, 1512, 1680, 1848, 2016, 2184, 2352, 2520, 2688, 2856])
ExRatio = np.array([2, 27, 43, 36, 39, 32, 27, 22, 13, 10, 14, 14, 4, 4, 7, 3, 3, 1])
ratio = 100* logistic_out[:,1]/(logistic_out[:,0]+logistic_out[:,1])
plt.plot(t,ratio)
plt.plot(time,ExRatio,'h')
xlabel('Position')
ylabel('Pollution')

New Cell

from scipy.stats import norm

# set up a function to return the log likelihood
def logistic_loglik(params,t,data,sig):
    return sum(norm.logpdf(logistic(data, t, params),sig))

# set standard deviations to be 10% of the population values
sig = ExRatio/10

# parameters for the MCMC
reps = 50000
npars = 3

# output matrix
par_out = np.ones(shape=(reps,npars))
# acceptance 
accept = np.zeros(shape=(reps,npars))
# proposal standard deviations. These have been pre-optimized.
propsigma = [0.05,20,5]

for i in range(1,reps):
    # make a copy of previous parameters
    par_out[i,] = par_out[i-1,]
    for j in range(npars):
        proposed = np.copy(par_out[i,:]) # we need to make a copy so that rejected moves don't affect the original matrix
        proposed[j] = proposed[j] + np.random.normal(0,propsigma[j])
        if (proposed[j]>0): # automatically reject moves if proposed parameter <=0 
            alpha = np.exp(logistic_loglik(proposed,time,ExRatio,sig)-logistic_loglik(par_out[i-1,],time,ExRatio,sig))
            u = np.random.rand()
            if (u < alpha):
                par_out[i,j] = proposed[j]
                accept[i,j] = 1

#print(sum(accept[range(101,reps),:])/(reps-100))


#plt.plot(par_out[:,0])
#plt.plot(par_out[range(101,reps),0])
#plt.plot(par_out[:,0],par_out[:,2])
plt.hist(par_out[range(101,reps),0],50)
print('\n')
a=np.mean(par_out[range(101,reps),0])

I think its mistaking my parameters for something else but that might be wrong. I am using Jupyter notebook

Answer 1

You cannot use S, R, A = x, if x is empty or has not enough (too much) values to unpack.

For what I see, you are trying to define S, R and A values using the variable x. It is possible in that way only if x is the len of 3. If you want to assign certain x values to specific S, R or A use loop, or if you want to do this that way you can use:

S, R, *A = x,

this way the variable S and R will have the first and second element of x, and variable A the rest. You can put * before any variable to make it take the excessive values you store in x.

Answer 2

The immediate syntax error

In your call

---> 28             alpha = np.exp(logistic_loglik(proposed,time,ExRatio,sig)-logistic_loglik(par_out[i-1,],time,ExRatio,sig))

to

      4 def logistic_loglik(params,t,data,sig):
----> 5     return sum(norm.logpdf(logistic(data, t, params),sig))

where finally the paramerers are used as defined in

      7 def logistic(x,t,params):
----> 8     S, R, A = x

the x that causes the error is the data of the previous call which is set to exTatio in the first call which is defined in your first block as an array of several handful of values. There might be something wrong with the logic that you use, as the structure of exRatio is not the one of 3 state variables.

---

The correct implementation of the concept of the log-likelihood sum

What you want is to compute the log-likelihood of the computed ratios at your sample points where the distribution for ExTime[k] is given as a normal distribution with mean ExRatio[k] and variance sig[k] which is set to ExRatio[k]/10. In your code you need to do exactly that, solve the ODE with the proposed initial values, compute the ratios and sum the logs of the pdf values:

# set up a function to return the log likelihood
def logistic_loglik(SRA0,ExTime,ExRatio,sig):
    # solve the ODE with the trial values `SRA0` and 
    # output the samples at the sample times `ExTime`
    logistic_out = odeint(logistic, SRA0, ExTime, args=(params,))
    # compute the ratios
    ratio = 100* logistic_out[:,1]/(logistic_out[:,0]+logistic_out[:,1])
    # return the summed log-likelihood
    return sum(norm.logpdf(ratio, ExRatio, sig))

Trying variants of propsigma leads to initially rapid convergence to qualitatively reasonable fits.

propsigma       i       S, R, A = par_out[i]

[0.05,20.,5.]  59 [ 2.14767909    0.18163897   5.45312544]
[20,0.5,5.]    39 [ 56.48959836   0.50890498   5.80229728]
[5.,2.,5.]     79 [ 67.26394337   0.15865463   6.0213663 ]