0

I have this example code that i am using in python with pymc3 package

import pymc3
import numpy as np


size = 100

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

X2 = np.random.randn(size) * 0.2

X3 = np.random.randn(size) * 2

alpha, sigma = 1, 1
beta = [1, 2.5]

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

basic_model = pymc3.Model()

with basic_model:

    # Priors for unknown model parameters
    alpha = pymc3.Normal('alpha', mu=0, sd=1)
    beta = pymc3.Normal('beta', mu=0, sd=1, shape=3)
    sigma = pymc3.HalfNormal('sigma', sd=1)
    # Expected value of outcome
    mu =  alpha + beta[0]*X1 + beta[1]*X2 + beta[2]*X3

    Y_obs = pymc3.Normal('Y_obs', mu=mu, sd=sigma, observed = Y)
    trace = pymc3.sample(100,chains = 1,  step = pymc3.NUTS())
    a = trace['alpha']
    b = trace['beta']
    print  (trace['alpha'])
    print (trace['beta'])
print (a[0] + b[0][0]*X1[0] + b[0][1]*X2[0] + b[0][2]*X3[0]) #mu

I expected that if mu = 0, only beta and alpha constants will vary? And I get an analysis from that, but mu result is never 0 or even close to 0? For example, if alpha is negative other constants or at least one should be positive so every the calculation(mu) should be 0 or close to 0?. If I try to set mu to = 0 I get an error SyntaxError: can't assign to literal

Noob Programmer
  • 698
  • 2
  • 6
  • 22
  • It is not clear what you are trying to achieve here... What is mu in your example? You are not conditioning on any value for mu and as such there is no likelihood term? So the values for a, b[:,0], b[:,1], b[:,2] are all draws from your standard normal prior and there is no other connection in the model... – nick Jun 06 '19 at 10:30
  • Sorry, I am new to this. I added some code. mu is that equation right? I wanna do an analysis to vary parameters . Like I said upside in the text. Beta and alpha constant need to vary based on X1, X2, X3 with the condition that mu is 0(or some other number). For example if X1,X2,X3 = 1, then if alpha = 3 -> b[0],b[1],b[2] = -1. So 0 = 3 + (-1*1) + (-1*1) + (-1*1) = 0 – Noob Programmer Jun 06 '19 at 10:45
  • So it is pretty uncommon to refer to "varying parameters". In the setup here we imagine that the parameters have a distribution (N(0,1)) and these, along with some features X1, X2, X3 will produce a response variable (mu) in your case. Are you sure pymc3 is the tool you want to use here? It seems like some sort of linear programming technique might be more helpful? – nick Jun 06 '19 at 10:55
  • In your example, you do realise that the solution is not unique right? b[0],b[1],b[2] = -2,-1,0 seems to be as plausible? – nick Jun 06 '19 at 10:57
  • Yeah sure, that was just an example for easier understanding of what I am trying to do. More options are plausible. I am trying to do some sort of sensitivity analysis fro Bayesian inference. So in my model I am gonna have around 1000 parameters and I wanna check which one changes if I change a specific one. So if I have paramaters [a1 a2.....an], how much and in what way will (will the value fall or rise) a2 change if I change value of a1. – Noob Programmer Jun 06 '19 at 11:03
  • 1
    Ah I see... So correct me if I am wrong: We have a model where (let's start calling mu 'y' so there is no confusion between 'mu' param of your Normal variable and mu, the observed data) y ~ a + b1*x1 + b2*x2 + b3*x3. Now you'd like to perform inference and say for example, conditioning on a > 3, what are likely values of b1, b2, b3? – nick Jun 06 '19 at 11:11
  • Note my use of a>3 there. You have assumed a is a continuous variable and so Pr(a=3) = 0. I.e, the probability that a = a specific value is 0. – nick Jun 06 '19 at 11:12
  • Yes for example :) and in what relations they are among themselves, so if I rise the value of one, the second one needs to lower for some amount. :) oh Thanks for your time! – Noob Programmer Jun 06 '19 at 11:14
  • Sure thing. Ok cool. I still just need to think through exactly what you are after and I'll try get around to posting a reply in a bit. – nick Jun 06 '19 at 11:23
  • Ah I see you included the `observed` data. This is much more helpful! – nick Jun 06 '19 at 11:24
  • 1
    Well, take it to the consideration that it could be wrong :( like I said I'm doing that for three/four days now, still a lot confused :( – Noob Programmer Jun 06 '19 at 11:26
  • Hi. I have taken a look at your question and I am afraid I still don't follow... I think you need to clarify what it is you are after. You are essentially running Bayesian linear regression in this example with 4 unknown parameters (a, b1, b2 and b3). pymc3 gives you posterior distributions over these parameters. I suppose you could put a strong prior over a parameter (say alpha?) to "change" it but again it is not clear to me why you would want to do that. Perhaps you are after the sensitivity of the model to changes in *data* in which case maybe you want to vary X and not the params? – nick Jun 06 '19 at 15:14
  • Yes exactly I am after the sensitivity of the model. How do I change my code to that? Since the pymc3 only calculates alpha, beta and sigma? – Noob Programmer Jun 07 '19 at 06:29
  • Can we go in a private chatroom? so we don't spamm here :) – Noob Programmer Jun 07 '19 at 08:31

0 Answers0