Hello,everyone.
I’m trying to carry out a sampling of the multiple Gauss mixed model，The distribution of variables is as follows(Mixed multivariate Gauss distribution)：
y_i=∑(k=1)^K{ω_kN(O,σ_k^2 H(ρ_k ))} ----multivariate Gauss distribution mixture_
ω_k=β_k∏(j=1)^(k-1){1-β_j } ----weights,stick-breaking process_
β_i iid Beta(1,M)
σ_k iid U(0,1.5) prior
ρ_k iid U(0,1) prior
I did the following the code but it is not working：
import numpy as np
import numpy as np
import pymc3 as pm
from theano import tensor as tt
K = 6 # the number of components
N = 200 # the number of observed individuals
n = 9 # each individual was observed 9 times
mu0 = np.linspace(0.0, 0.0, num=n)
# simulate observation data
C1 = np.zeros((n, n))
for i in range(0, n):
for j in range(0, n):
if i == j:
C1[i, j] = 10.0
else:
C1[i, j] = 7.0
dataSet1 = np.random.multivariate_normal(mu0, C1, size=N)
# observation time(all individuals are observed at the same times)
time_obseved = [-5.0, -4.5, -3.0, -2.5, -1.0, 0.0, 1.5, 2.5, 3.0, 4.0, 5.5]
time = np.array(time_obseved)
#I missed it before, and I'm mending it now
def H(rho):
H=np.zeros((n,n))
for i in range(0,n):
for j in range(0,n):
H[i,j]=np.power(rho,np.abs(time[j]-time[i]))
return H
# stick-breaking process
def stick_breaking(beta):
portion_remaining = tt.concatenate(
[[1], tt.extra_ops.cumprod(1 - beta)[:-1]])
return beta * portion_remaining
with pm.Model() as model:
M = pm.Gamma('M', 1., 1.)
sigma_w = pm.Uniform('sigma_w', 0.0, 1.5, shape=K)
rho = pm.Uniform('rho', 0.0, 0.1, shape=K)
beta = pm.Beta('beta', 1., M, shape=K)
w = pm.Deterministic('w', stick_breaking(beta))
obs = pm.Mixture('obs', w, pm.MvNormal.dist(
mu=mu0, cov=sigma_w**2 * H(rho)), observed=dataSet1)
I get the following error:
runfile('F:/spyder/Mixed multivariate Gauss distribution.py')
WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\Administrator\Anaconda3\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory
warnings.warn("DeprecationWarning: there is no c++ compiler."
WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string.
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
C:\Users\Administrator\Anaconda3\lib\site-packages\h5py\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
from ._conv import register_converters as _register_converters
Traceback (most recent call last):
File "<ipython-input-1-a62ff05957e7>", line 1, in <module>
runfile('F:/spyder/Mixed multivariate Gauss distribution.py')
File "C:\Users\Administrator\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 705, in runfile
execfile(filename, namespace)
File "C:\Users\Administrator\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 102, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)
File "F:/spyder/Mixed multivariate Gauss distribution.py", line 75, in <module>
obs=pm.Mixture('obs',w,pm.MvNormal.dist(mu=mu0,cov=sigma_w**2*H(rho)),observed=dataSet1)
File "F:/spyder/Mixed multivariate Gauss distribution.py", line 54, in H
H[i,j]=np.power(rho,np.abs(time[j]-time[i]))
ValueError: setting an array element with a sequence.
How could I make this work?
