pymc.Normal object initialised with single values, returns an array of 184 values.

[mar-ses]

Hello, I've come across a problem with some MCMC modelling I'm doing. I've got 184 individual planets, each one with a true radius R, mass M, and insolation flux F. I have a forward model with calculates R, given M, F and the model hyperparameters.

Recently I've changed something in my model and it's stopped working. I cannot debug it or figure out why. Here is the relevant section of code

R = np.empty(N, dtype=object)

for i in range(0,N):
    print("mu data type: " + str(np.result_type(R_0 + A*(abs(F[i])**gamma) + B_array[i]*M[i])) )
    print("tau data type: " + str(np.result_type(tauR)) )

    R[i] = pymc.Normal('R_%i' %i, mu = (R_0 + A*(abs(F[i])**gamma) + B_array[i]*M[i]), tau=tauR)

#               DATA PARAMETERS

print("R length = " + str(len(R)) )
print("R_1 length = " + str(len(R[1].value)) )
print("R_1 value = " + str(R[1].value) )

R_obs = pymc.Normal('Observed radius', mu = R, tau = 1/(R_stdev**2), value = R_obs_array, observed = True)

The print statements are there for debugging.

So basically, R is an array of N Stochastic Objects (pymc.Normal). Its mu is calculated from the statement:

R_0 + A*(abs(F[i])**gamma) + B_array[i]*M[i]

I have checked (in the print statement): this statement is definitely a float64, not an array. Its tau is tauR, again 100% sure that it's a single float.

Now we get to the next section; len(R) = N as expected.

However, suddently, len(R[1].value) = 184 as well. I can't see how that would come about, all the arguments for creating R[1] as a pymc.Normal object are floats. I can't see R[i].value return 184 values, instead of 1. Can anyone help me here?

My model was working previously, with a similar situation except without the parameters M and B_array.

[mar-ses]

Here is an extract of the previous working model, it's got a very similar situation, in the above I tried to streamline it and make every individual part work because it was giving me the error.

@pymc.deterministic(plot=False)
def mu(F=F, A=A, g=gamma, R_0=R_0 ):
    return R_0 + A*(abs(F)**g)

@pymc.deterministic(plot=False)
def tau( log_var = log_var ):
    return 1/math.exp(log_var)

R = np.empty(N, dtype=object)

for i in range(0,N):
    R[i] = pymc.Normal('R_%i' %i, mu[i], tau)

#               DATA PARAMETERS

print("R length = " + str(len(R)) )
print("R_obs length = " + str(len(R_obs_array)) )
print("R_std length = " + str(len(R_stdev)) )
print("mu length = " + str(len(mu.value)) + ' ' + str(len(mu[1].value)) )

R_obs = pymc.Normal('Observed radius', mu = R, tau = 1/(R_stdev**2), value = R_obs_array, observed = True)

The difference here is: there is no masses M, and no B_array. So it calculates an array of mu's (the mean radius for each planet) from the array F and the hyperparameters. Again, R is an array of N objects, with R[i] = pymc.Normal(..., mu[i], tau[i] ).

As far as I can tell it's the exact same situation, except here I get what I want, which is that R[i] is not an array of 184 floats or whatever, unlike in the OP.

[fonnesbeck]

Does this mean you solved the issue?