问题描述
我在尝试对模型进行编码的代码中遇到问题，出现了以下错误，作为一个相对新手，我不确定如何解决.
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)
我尝试编写的模型是一个MCMC，以使某些ODE适应某些数据，我在下面为上下文添加了代码.
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')
新单元格
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])
我认为它误将我的参数用于其他方面，但这可能是错误的.我正在使用Jupyter笔记本
I think its mistaking my parameters for something else but that might be wrong.I am using Jupyter notebook
推荐答案
立即语法错误
在通话中
---> 28 alpha = np.exp(logistic_loglik(proposed,time,ExRatio,sig)-logistic_loglik(par_out[i-1,],time,ExRatio,sig))
到
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
导致错误的
x是上一个调用的data，该调用在第一个调用中被设置为exTatio，在您的第一个块中将其定义为几个值的数组.您使用的逻辑可能有问题，因为exRatio的结构不是3个状态变量之一.
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.
您想要的是在样本点上计算计算比率的对数似然，其中ExTime[k]的分布为正态分布，均值ExRatio[k]且方差sig[k]设置为.在您的代码中，您需要执行此操作，使用建议的初始值求解ODE，计算比率并将pdf值的对数求和:
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))
尝试使用propsigma的变体会导致最初的快速收敛到定性的合理拟合.
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 ]
这篇关于ValueError:太多值无法解包(预期3)?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！