Does having had a visa for a country mean I used to be a citizen/national of that country?

Why are < or > required to use /dev/tcp

Methodology: Writing unit tests for another developer

Count All Possible Unique Combinations of Letters in a Word

How does the spell Remove Curse interact with a Sword of Vengeance?

Why tighten down in a criss-cross pattern?

Is "Busen" just the area between the breasts?

Can a rogue use Sneak Attack in a Darkness spell cast by another player?

How to find the last non zero element in every column throughout dataframe?

If plants "alternate generations" between sporophytes and gametophytes, why don't we say the same of humans?

Greeting with "Ho"

Why is prior to creation called holy?

Do I need a shock-proof watch for cycling?

When two first person POV characters meet

What does the hyphen "-" mean in "tar xzf -"?

How many people are necessary to maintain modern civilisation?

Helping ease my back pain when I'm studying 13 hours everyday, even weekends

What exactly is the 'online' in OLAP and OLTP?

Prime sieve in Python

Concurrent normals conjecture

Does a vocal melody have any rhythmic responsibility to the underlying arrangement in pop music?

What does "play with your toy’s toys" mean?

Are all instances of trolls turning to stone ultimately references back to Tolkien?

How does a blind passenger not die, if driver becomes unconscious

PyMC3 Normal with variance per column

Custom likelihood in pymc3Using a complex likelihood in PyMC3PyMC3 - Index 2-dimensional data while fitting hierarchical autoregressive modelHow to write a custom Deterministic or Stochastic in pymc3 with theano.op?how to fit a method belonging to an instance with pymc3?sampling multivariate uniform in PyMC3Multiple levels in hierarchical linear regression using PYMC3Fusing Multiple Measurements in pyMC3Comparing models fitting multivariate dataSlow, underlying trend with a random walk or Gaussian process

.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty height:90px;width:728px;box-sizing:border-box;

1

I am trying to define a pymc3.Normal variable with the following as mu:

import numpy as np
import pymc3 as pm

mx = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)

The problem is that all random variables in p_ get the same std (epsilon). I would like the first row to use epsilon1, the second row epsilon2 etc.

How Can I do that?

pymc3 pymc

share|improve this question

edited Mar 25 at 8:49

Eran

asked Mar 25 at 8:23

EranEran

182110

add a comment | 

1

I am trying to define a pymc3.Normal variable with the following as mu:

import numpy as np
import pymc3 as pm

mx = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)

The problem is that all random variables in p_ get the same std (epsilon). I would like the first row to use epsilon1, the second row epsilon2 etc.

How Can I do that?

pymc3 pymc

share|improve this question

edited Mar 25 at 8:49

Eran

asked Mar 25 at 8:23

EranEran

182110

add a comment | 

I am trying to define a pymc3.Normal variable with the following as mu:

import numpy as np
import pymc3 as pm

mx = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)

The problem is that all random variables in p_ get the same std (epsilon). I would like the first row to use epsilon1, the second row epsilon2 etc.

How Can I do that?

pymc3 pymc

share|improve this question

edited Mar 25 at 8:49

Eran

asked Mar 25 at 8:23

EranEran

182110

import numpy as np
import pymc3 as pm

mx = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)

share|improve this question

edited Mar 25 at 8:49

Eran

asked Mar 25 at 8:23

EranEran

182110

share|improve this question

edited Mar 25 at 8:49

Eran

asked Mar 25 at 8:23

EranEran

182110

asked Mar 25 at 8:23

EranEran

182110

asked Mar 25 at 8:23

EranEran

182110

1 Answer
1

active

oldest

votes

0

One can pass an argument for the shape parameter to achieve this. To demonstrate this, let's make some fake data to pass as observed, where we use fixed values for epsilon that we can compare against the inferred ones.

Example Model

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
 # use a (1,4) matrix to pool variance by columns
 epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

 p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

 trace = pm.sample(random_seed=2019)

This samples well, and gives the following summary

which clearly bound the true values of the standard deviations within the HPDs.

share|improve this answer

answered Apr 1 at 20:44

mervmerv

27.3k676115

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  The sd is sd=epsilon which is of shape (1, mu.shape[1]). This does the trick?
  
  – Eran
  Apr 6 at 16:18
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Eran yes, exactly
  
  – merv
  Apr 6 at 16:21
  

  
  
  
  

add a comment | 

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

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f55333706%2fpymc3-normal-with-variance-per-column%23new-answer', 'question_page');

);

Post as a guest

Name

Email

Required, but never shown

0

One can pass an argument for the shape parameter to achieve this. To demonstrate this, let's make some fake data to pass as observed, where we use fixed values for epsilon that we can compare against the inferred ones.

Example Model

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
 # use a (1,4) matrix to pool variance by columns
 epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

 p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

 trace = pm.sample(random_seed=2019)

This samples well, and gives the following summary

which clearly bound the true values of the standard deviations within the HPDs.

share|improve this answer

answered Apr 1 at 20:44

mervmerv

27.3k676115

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  The sd is sd=epsilon which is of shape (1, mu.shape[1]). This does the trick?
  
  – Eran
  Apr 6 at 16:18
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Eran yes, exactly
  
  – merv
  Apr 6 at 16:21
  

  
  
  
  

add a comment | 

0

One can pass an argument for the shape parameter to achieve this. To demonstrate this, let's make some fake data to pass as observed, where we use fixed values for epsilon that we can compare against the inferred ones.

Example Model

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
 # use a (1,4) matrix to pool variance by columns
 epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

 p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

 trace = pm.sample(random_seed=2019)

This samples well, and gives the following summary

which clearly bound the true values of the standard deviations within the HPDs.

share|improve this answer

answered Apr 1 at 20:44

mervmerv

27.3k676115

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  The sd is sd=epsilon which is of shape (1, mu.shape[1]). This does the trick?
  
  – Eran
  Apr 6 at 16:18
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @Eran yes, exactly
  
  – merv
  Apr 6 at 16:21
  

  
  
  
  

add a comment | 

One can pass an argument for the shape parameter to achieve this. To demonstrate this, let's make some fake data to pass as observed, where we use fixed values for epsilon that we can compare against the inferred ones.

Example Model

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
 # use a (1,4) matrix to pool variance by columns
 epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

 p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

 trace = pm.sample(random_seed=2019)

This samples well, and gives the following summary

which clearly bound the true values of the standard deviations within the HPDs.

share|improve this answer

answered Apr 1 at 20:44

mervmerv

27.3k676115

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ], 
 [0.25 , 0.333, 0.25 , 0. ],
 [0.25 , 0.167, 0. , 0. ],
 [0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
 # use a (1,4) matrix to pool variance by columns
 epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

 p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

 trace = pm.sample(random_seed=2019)

share|improve this answer

answered Apr 1 at 20:44

mervmerv

27.3k676115

answered Apr 1 at 20:44

mervmerv

27.3k676115

answered Apr 1 at 20:44

mervmerv

27.3k676115