PyMC3 passing stochastic covariance matrix to pm.MvNormal()










1















I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.



import numpy as np
import pymc3 as pm

n = 10000;
np.random.seed(0)
X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

with pm.Model() as model:
# PRIORS
mu = [pm.Uniform('mux', lower=-1, upper=1),
pm.Uniform('muy', lower=-1, upper=1)]
cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
[0, pm.Uniform('a22', lower=0.1, upper=2)]])

# LIKELIHOOD
likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

with model:
trace = pm.sample(draws=1000, chains=2, tune=1000)


while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.



Where am I going wrong?










share|improve this question




























    1















    I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.



    import numpy as np
    import pymc3 as pm

    n = 10000;
    np.random.seed(0)
    X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

    with pm.Model() as model:
    # PRIORS
    mu = [pm.Uniform('mux', lower=-1, upper=1),
    pm.Uniform('muy', lower=-1, upper=1)]
    cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
    [0, pm.Uniform('a22', lower=0.1, upper=2)]])

    # LIKELIHOOD
    likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

    with model:
    trace = pm.sample(draws=1000, chains=2, tune=1000)


    while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.



    Where am I going wrong?










    share|improve this question


























      1












      1








      1








      I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.



      import numpy as np
      import pymc3 as pm

      n = 10000;
      np.random.seed(0)
      X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

      with pm.Model() as model:
      # PRIORS
      mu = [pm.Uniform('mux', lower=-1, upper=1),
      pm.Uniform('muy', lower=-1, upper=1)]
      cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
      [0, pm.Uniform('a22', lower=0.1, upper=2)]])

      # LIKELIHOOD
      likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

      with model:
      trace = pm.sample(draws=1000, chains=2, tune=1000)


      while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.



      Where am I going wrong?










      share|improve this question
















      I've tried to fit a simple 2D gaussian model to observed data by using PyMC3.



      import numpy as np
      import pymc3 as pm

      n = 10000;
      np.random.seed(0)
      X = np.random.multivariate_normal([0,0], [[1,0],[0,1]], n);

      with pm.Model() as model:
      # PRIORS
      mu = [pm.Uniform('mux', lower=-1, upper=1),
      pm.Uniform('muy', lower=-1, upper=1)]
      cov = np.array([[pm.Uniform('a11', lower=0.1, upper=2), 0],
      [0, pm.Uniform('a22', lower=0.1, upper=2)]])

      # LIKELIHOOD
      likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)

      with model:
      trace = pm.sample(draws=1000, chains=2, tune=1000)


      while I can do this in 1D by passing the sd to pm.Normal, I have some trouble in passing the covariance matrix to pm.MvNormal.



      Where am I going wrong?







      python theano bayesian normal-distribution pymc3






      share|improve this question















      share|improve this question













      share|improve this question




      share|improve this question








      edited Nov 14 '18 at 22:42









      merv

      25.2k673109




      25.2k673109










      asked Nov 9 '18 at 9:05









      FabioFabio

      1335




      1335






















          1 Answer
          1






          active

          oldest

          votes


















          1














          PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.



          Original Model



          Keeping with the intent of your code, a working version would go something like



          import numpy as np
          import pymc3 as pm
          import theano.tensor as tt

          N = 10000
          np.random.seed(0)
          X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

          with pm.Model() as model:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # diagonal values on covariance matrix
          a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

          # convert vector to a 2x2 matrix with `a` on the diagonal
          cov = tt.diag(a)

          likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)


          Alternative Model



          I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.



          Here's a simple application of that example in this context:



          with pm.Model() as model_lkj:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # LKJ prior for covariance matrix (see example)
          packed_L = pm.LKJCholeskyCov('packed_L', n=2,
          eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
          # convert to (2,2)
          L = pm.expand_packed_triangular(2, packed_L)

          likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)





          share|improve this answer

























          • merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

            – Fabio
            Nov 15 '18 at 18:32










          Your Answer






          StackExchange.ifUsing("editor", function ()
          StackExchange.using("externalEditor", function ()
          StackExchange.using("snippets", function ()
          StackExchange.snippets.init();
          );
          );
          , "code-snippets");

          StackExchange.ready(function()
          var channelOptions =
          tags: "".split(" "),
          id: "1"
          ;
          initTagRenderer("".split(" "), "".split(" "), channelOptions);

          StackExchange.using("externalEditor", function()
          // Have to fire editor after snippets, if snippets enabled
          if (StackExchange.settings.snippets.snippetsEnabled)
          StackExchange.using("snippets", function()
          createEditor();
          );

          else
          createEditor();

          );

          function createEditor()
          StackExchange.prepareEditor(
          heartbeatType: 'answer',
          autoActivateHeartbeat: false,
          convertImagesToLinks: true,
          noModals: true,
          showLowRepImageUploadWarning: true,
          reputationToPostImages: 10,
          bindNavPrevention: true,
          postfix: "",
          imageUploader:
          brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
          contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
          allowUrls: true
          ,
          onDemand: true,
          discardSelector: ".discard-answer"
          ,immediatelyShowMarkdownHelp:true
          );



          );













          draft saved

          draft discarded


















          StackExchange.ready(
          function ()
          StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53222664%2fpymc3-passing-stochastic-covariance-matrix-to-pm-mvnormal%23new-answer', 'question_page');

          );

          Post as a guest















          Required, but never shown

























          1 Answer
          1






          active

          oldest

          votes








          1 Answer
          1






          active

          oldest

          votes









          active

          oldest

          votes






          active

          oldest

          votes









          1














          PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.



          Original Model



          Keeping with the intent of your code, a working version would go something like



          import numpy as np
          import pymc3 as pm
          import theano.tensor as tt

          N = 10000
          np.random.seed(0)
          X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

          with pm.Model() as model:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # diagonal values on covariance matrix
          a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

          # convert vector to a 2x2 matrix with `a` on the diagonal
          cov = tt.diag(a)

          likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)


          Alternative Model



          I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.



          Here's a simple application of that example in this context:



          with pm.Model() as model_lkj:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # LKJ prior for covariance matrix (see example)
          packed_L = pm.LKJCholeskyCov('packed_L', n=2,
          eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
          # convert to (2,2)
          L = pm.expand_packed_triangular(2, packed_L)

          likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)





          share|improve this answer

























          • merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

            – Fabio
            Nov 15 '18 at 18:32















          1














          PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.



          Original Model



          Keeping with the intent of your code, a working version would go something like



          import numpy as np
          import pymc3 as pm
          import theano.tensor as tt

          N = 10000
          np.random.seed(0)
          X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

          with pm.Model() as model:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # diagonal values on covariance matrix
          a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

          # convert vector to a 2x2 matrix with `a` on the diagonal
          cov = tt.diag(a)

          likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)


          Alternative Model



          I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.



          Here's a simple application of that example in this context:



          with pm.Model() as model_lkj:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # LKJ prior for covariance matrix (see example)
          packed_L = pm.LKJCholeskyCov('packed_L', n=2,
          eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
          # convert to (2,2)
          L = pm.expand_packed_triangular(2, packed_L)

          likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)





          share|improve this answer

























          • merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

            – Fabio
            Nov 15 '18 at 18:32













          1












          1








          1







          PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.



          Original Model



          Keeping with the intent of your code, a working version would go something like



          import numpy as np
          import pymc3 as pm
          import theano.tensor as tt

          N = 10000
          np.random.seed(0)
          X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

          with pm.Model() as model:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # diagonal values on covariance matrix
          a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

          # convert vector to a 2x2 matrix with `a` on the diagonal
          cov = tt.diag(a)

          likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)


          Alternative Model



          I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.



          Here's a simple application of that example in this context:



          with pm.Model() as model_lkj:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # LKJ prior for covariance matrix (see example)
          packed_L = pm.LKJCholeskyCov('packed_L', n=2,
          eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
          # convert to (2,2)
          L = pm.expand_packed_triangular(2, packed_L)

          likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)





          share|improve this answer















          PyMC3 distribution objects are not simple numeric objects or numpy arrays. Instead, they are nodes in a theano computation graph and often require operations from either pymc3.math or theano.tensor to manipulate them. Moreover, placing PyMC3 objects in numpy arrays is unnecessary since they already are multidimensional.



          Original Model



          Keeping with the intent of your code, a working version would go something like



          import numpy as np
          import pymc3 as pm
          import theano.tensor as tt

          N = 10000
          np.random.seed(0)
          X = np.random.multivariate_normal(np.zeros(2), np.eye(2), size=N)

          with pm.Model() as model:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # diagonal values on covariance matrix
          a = pm.Uniform('a', lower=0.1, upper=2, shape=2)

          # convert vector to a 2x2 matrix with `a` on the diagonal
          cov = tt.diag(a)

          likelihood = pm.MvNormal('likelihood', mu=mu, cov=cov, observed=X)


          Alternative Model



          I assume the example you provided is just a toy to communicate the problem. But just in case, I'll mention that the main advantage of using a multivariate normal (modeling covariance between parameters) is lost when restricting the covariance matrix to be diagonal. Furthermore, the theory of priors for covariance matrices is well-developed, so it's worth one's time considering existing solutions. In particular, there is a PyMC3 example using the LKJ prior for covariance matrices.



          Here's a simple application of that example in this context:



          with pm.Model() as model_lkj:
          # use `shape` argument to define tensor dimensions
          mu = pm.Uniform('mu', lower=-1, upper=1, shape=2)

          # LKJ prior for covariance matrix (see example)
          packed_L = pm.LKJCholeskyCov('packed_L', n=2,
          eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
          # convert to (2,2)
          L = pm.expand_packed_triangular(2, packed_L)

          likelihood = pm.MvNormal('likelihood', mu=mu, chol=L, observed=X)






          share|improve this answer














          share|improve this answer



          share|improve this answer








          edited Nov 14 '18 at 21:59

























          answered Nov 14 '18 at 20:40









          mervmerv

          25.2k673109




          25.2k673109












          • merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

            – Fabio
            Nov 15 '18 at 18:32

















          • merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

            – Fabio
            Nov 15 '18 at 18:32
















          merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

          – Fabio
          Nov 15 '18 at 18:32





          merv, thanks for the answer. Mine example was a toy model and I could use other ways as you suggested. In general I have difficulties in understanding how to include in the model custom operators. I have a topic similar on discourse.pymc.io/t/… where I have problem in including a "physical law" for the stochastic mean.

          – Fabio
          Nov 15 '18 at 18:32



















          draft saved

          draft discarded
















































          Thanks for contributing an answer to Stack Overflow!


          • Please be sure to answer the question. Provide details and share your research!

          But avoid


          • Asking for help, clarification, or responding to other answers.

          • Making statements based on opinion; back them up with references or personal experience.

          To learn more, see our tips on writing great answers.




          draft saved


          draft discarded














          StackExchange.ready(
          function ()
          StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53222664%2fpymc3-passing-stochastic-covariance-matrix-to-pm-mvnormal%23new-answer', 'question_page');

          );

          Post as a guest















          Required, but never shown





















































          Required, but never shown














          Required, but never shown












          Required, but never shown







          Required, but never shown

































          Required, but never shown














          Required, but never shown












          Required, but never shown







          Required, but never shown







          Popular posts from this blog

          Top Tejano songwriter Luis Silva dead of heart attack at 64

          政党

          天津地下鉄3号線