Any suggestion is appreciated.
I run the code listed below and obtain error the message “Size length is incompatible with batched dimensions of parameter 0 alpha: len(size) = 1, len(batched dims alpha) = 2. Size length must be 0 or >= 2”. My code is listed below. If I changed the link function
mu = pm.Deterministic("mu", pm.math.invlogit(a*np.sqrt(x)+b))
to mu = pm.Deterministic("mu", pm.math.invlogit(a[0] + a[1] * x[:, 1]))
It works fine. However it is not what I want. Thank you very much for your any thoughts.
> RANDOM_SEED = 8927
> rng = np.random.default_rng(RANDOM_SEED)
>
> data = pd.read_excel("BN_samples1.xlsx", usecols=['Effective water level above ground [m]', 'Building damage ratio'])
> sns.pairplot(data=data, kind="scatter");
> data["Intercept"] = np.ones(len(data))
> labels = ["Intercept", "Effective water level above ground [m]"]
> x = data[labels].to_numpy()
>
> rloss = data['Building damage ratio'].to_numpy()
> var = data['Effective water level above ground [m]'].var()
>
> #train-test split
> test_indices = range(len(x) - 10, len(x))
> train_indices = range(len(x) - 10)
> x_train, x_test = x[train_indices], x[test_indices]
> rloss_train, rloss_test = rloss[train_indices], rloss[test_indices]
>
> #Define and fit the model
> coords = {"coeffs": labels}
>
> with pm.Model(coords=coords) as model:
> # Data containers
> x = pm.MutableData("x", x_train)
> rloss = pm.MutableData("rloss", rloss_train)
>
> # Priors
> a = pm.Normal("a", mu=0, sigma=1, dims="coeffs")
> b = pm.Normal("b", mu=0, sigma=1, dims="coeffs")
>
> # Link function
> mu = pm.Deterministic("mu", pm.math.invlogit(a*np.sqrt(x)+b))
> phi = pm.Deterministic("phi", pm.math.exp(mu*(1-mu)/var-1))
>
> alpha = pm.Deterministic("alpha", mu * phi)
> beta = pm.Deterministic("beta", (1 - mu) * phi)
>
>
> def logp_beta(obsLoss, alpha, beta):
> return pm.logp(pm.Beta.dist(alpha=alpha, beta=beta), obsLoss)
>
> def random(alpha, beta, rng= None, size=None) :
> return scipy.stats.beta(alpha, beta).rvs(size)
>
> rloss_custom = pm.CustomDist('rloss_custom', alpha, beta, logp=logp_beta, random=random, observed=rloss, shape=x.shape[0])
>
> with model:
> step = pm.Metropolis()
> idata = pm.sample(4800, step=step, chains=4, cores=1)