How to get a value in an array based on the sampled value

Questions
1

i have following code:

# -*- coding: utf-8 -*-
"""
Created on Thu Mar  4 21:38:45 2021

@author: hysplit
"""
import arviz as az
import matplotlib.pyplot as plt
import numpy as np

%config InlineBackend.figure_format = 'retina'
# Initialize random number generator
RANDOM_SEED = 8927
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")

# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# Size of dataset
size = 100

# Predictor variable
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2

# Simulate outcome variable
Y = alpha + beta[0] * X1 + beta[1] * X2 + np.random.randn(size) * sigma

fig, axes = plt.subplots(1, 2, sharex=True, figsize=(10, 4))
axes[0].scatter(X1, Y, alpha=0.6)
axes[1].scatter(X2, Y, alpha=0.6)
axes[0].set_ylabel("Y")
axes[0].set_xlabel("X1")
axes[1].set_xlabel("X2");

import pymc3 as pm

print(f"Running on PyMC3 v{pm.__version__}")

basic_model = pm.Model()

m_test=[1,1,1,1,3,3,4,3]


with basic_model:

    # Priors for unknown model parameters
    alpha = pm.Uniform("alpha",  lower=0, upper=5)   
    beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Expected value of outcome
    
    #mu = alpha + beta[0] * X1 + beta[1] * X2
    mu = m_test[abs(int(alpha))] + beta[0] * X1 + beta[1] * X2
    #mu =  beta_1[0] * X1 + beta_1[1] * X2
    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)

map_estimate = pm.find_MAP(model=basic_model)
map_estimate

I want get a value from the array m_test based on the sampled alpha
but the code did not work,for it is wrong.
can anyone give me some suggestion?
thanks a lot.

1 Like
2

maybe what i am asking is not clear.

what i what to do is pass the value of the sampled parameter to an array.
what should i do?
thanks

1 Like
3

pm.find_MAP() will return you a single value for each parameter (the most likely a posteriori). You probably want to do pm.sample() which will return a bunch of sampled values from the posterior of each parameter.

1 Like
4

the line: m_test[abs(int(alpha))] is what i want to get.
but this is wrong。

X1 = np.random.randn(size)
basic_model = pm.Model()
m_test=[1,1,1,1,3,3,4,3]
with basic_model:
 # Priors for unknown model parameters
    alpha = pm.Uniform("alpha",  lower=0, upper=5,transform=None)   
  # Expected value of outcome
    mu = m_test[abs(int(alpha))]*X1
 # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal("Y_obs", mu=mu, observed=Y)
    trace = pm.sample(1000)
1 Like
5

Oh sorry I missed what your problem was. So you are not succeeding in indexing m_test with a random variable. What error are you seeing?

1 Like
6

@aspen I took the liberty of editing the posts to format the code in monospace and no markdown formatting. You can format text as code inline with one backtick: The word `code` will be rendered in markdown as plain monospace text, like so: “The word code will be rendered in markdown as plain monospace text”.

For code block, you have to use three backticks, the following markdown for example:

Here is my example code:

```
import numpy as np
x = np.linspace(0, 1)
```

gets rendered as:

Here is my example code:

import numpy as np
x = np.linspace(0, 1)

Welcome to PyMC Discourse!

7

when i use the DiscreteUniform distribution

    alpha = pm.DiscreteUniform("alpha",  lower=0, upper=5,transform=None)  
    mu = m_test[alpha]*X1

I get the error:

TypeError: list indices must be integers or slices, not FreeRV

When i use the uniform distribution

alpha = pm.Uniform("alpha",  lower=0, upper=5,transform=None)  
 mu = m_test[abs(int(alpha))]*X1

I get the error

int() argument must be a string, a bytes-like object or a number, not 'FreeRV'

So the problem is how I can index m_test based on the variable of FreeRV format.
thanks

8

ok
Now i know.
thanks for the editing and instruction.

9

You can use pm.Data, which takes care of converting your array and allows for indexing:

In [2]: X1 = np.random.randn(10) 
   ...: m_test=np.array([1,1,1,1,3,3,4,3]) 
   ...: Y = np.zeros(10) 
   ...: with pm.Model() as basic_model: 
   ...:     m_test_model = pm.Data('m_test_model', m_test) 
   ...:     alpha = pm.DiscreteUniform("alpha",  lower=0, upper=9)    
   ...:     mu = m_test_model[alpha]*X1 
   ...:     Y_obs = pm.Normal("Y_obs", mu=mu, observed=Y)

Alternatively to convert to int from a Uniform this might work:

alpha = pm.Uniform('alpha',  0, 9)
mu = m_test_model[pm.theanof.intX(alpha)] * X1

I can get the MAP, but normal sampling raises errors on my side…

10

It worked.
thanks for your help.