I’m trying to figure out how to fit a beta distribution to data, and to capture the parameters’ correlation.

(My ultimate goal is successive updating of the posteriors with new data. And at each update I want to be able to specify a correlated prior based on the correlation of the previous iteration’s posterior.)

A) What I * have* been able to do (with two models):

- Generate fake data from a beta distribution
- Fit a model to that data
- Extract (a, b) samples from the fitted model
- Fit a second model to the (a, b) samples that captures the samples’ correlation

B) What I’d * like* to be able to do (with a single model):

- Generate fake data from a beta distribution
- Fit a (single) model to that data that captures the parameters’ correlation and that allows me to specify a prior on that correlation.

In particular, my two-model (A) code (indirectly) captures the parameters’ strong positive correlation of 0.83 (in the second model within (A)):

```
mean sd
aMu 4.3672460 0.014254404
bMu 2.2089633 0.006574057
Rho[1,2] 0.8313109 0.015791045
```

… whereas my single-model (B) code doesn’t:

```
mean sd
a 4.3207410 0.2802322
b 2.1848314 0.1363592
Rho[1,2] -0.1489613 0.5907234
```

Presumably (?) this is because in (B) I’m fitting a single (a, b) pair to all the data (albeit with uncertainty about a and b), for which parameter correlation doesn’t even make sense.

I’m not sure if my single-model “I’d like to” goal (B) makes sense, so I’m hoping for some feedback/guidance.

Here’s my (A) code:

```
library(rstan)
options(mc.cores = parallel::detectCores())
## A.1) Generate fake data from a beta distribution
N = 400
beta_a = 4
beta_b = 2
betaSamps = rbeta( n=N, shape1=beta_a, shape2=beta_b )
dataListA1 = list(samp=betaSamps, N=N)
mA1 = "
data{
int N;
vector[N] samp;
}
parameters{
real<lower = 0> a;
real<lower = 0> b;
}
model{
samp ~ beta(a, b);
a ~ normal(2, 2);
b ~ normal(2, 2);
}
"
## A.2) Fit a model to that data
mA1Fit = stan( model_code=mA1, data=dataListA1 )
## A.3) Extract (a, b) samples from the fitted model
modelSamps = extract(mA1Fit)
## A.4) Fit a second model to the (a, b) samples that captures the samples'
# correlation
# No need for the full 4000 samples, so use just N
dataListA2 = list( aSamp=modelSamps$a[1:N], bSamp=modelSamps$b[1:N],
N=N )
mA2 = "
data{
int N;
vector[N] bSamp;
vector[N] aSamp;
}
parameters{
real aMu;
real bMu;
real<lower=0> bSigma;
real<lower=0> aSigma;
corr_matrix[2] Rho;
}
model{
bMu ~ normal( 3 , 2 );
aMu ~ normal( 4 , 2 );
bSigma ~ normal( 1, 1);
aSigma ~ normal( 1, 1);
Rho ~ lkj_corr( 1 );
{
vector[2] AB[N];
vector[2] MU;
MU = [ aMu , bMu ]';
vector[2] Sigma;
Sigma = [ aSigma , bSigma ]';
for ( j in 1:N ) AB[j] = [ aSamp[j] , bSamp[j] ]';
AB ~ multi_normal( MU , quad_form(Rho , diag_matrix(Sigma) ) );
}
}
"
fitMA2 = stan( model_code=mA2, data=dataListA2 )
fitMA2Summary = summary( fitMA2, pars = c("aMu", "bMu", "Rho[1,2]") )$summary
fitMA2Summary = fitMA2Summary[, c("mean", "sd")]
print(fitMA2Summary)
```

And here’s my unsuccessful single-model (B) code:

```
## B.1) Generate fake data from a beta distribution
# (Use data from step A1)
## B.2) Fit a model to that data that captures the parameters’ correlation and
# that allows me to specify a prior on that correlation.
mB = "
data{
int N;
vector[N] samp;
}
parameters{
real<lower = 0> a;
real<lower = 0> b;
corr_matrix[2] Rho;
}
model{
samp ~ beta(a, b);
Rho ~ lkj_corr( 1 );
{
vector[2] AB;
vector[2] MU;
MU = [ 4 , 3 ]';
vector[2] Sigma;
Sigma = [ 1 , 2 ]';
AB = [ a , b ]';
AB ~ multi_normal( MU , quad_form(Rho , diag_matrix(Sigma) ) );
}
}
"
fitMB = stan( model_code=mB, data=dataListA1 )
fitMBSummary = summary( fitMB, pars=c("a", "b", "Rho[1,2]") )$summary
fitMBSummary = fitMBSummary[,c("mean", "sd")]
print(fitMBSummary)
```