# Sufficient Statistics for Normal Case

**URL:** <https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171>\
**Category:** Modeling\
**Tags:** fitting-issues, specification, hierarchical-model\
**Created:** [June 24, 2021, 12:25pm UTC](https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171 "2021-06-24T12:25:27Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![eric11230](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric11230/32/11967_2.png) [@eric11230](https://discourse.mc-stan.org/u/eric11230)\
**Post date:** [June 24, 2021, 12:25pm UTC](https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171/1 "2021-06-24T12:25:27Z")

</div>

Hello,  
I do not think I have seen this solved but new to Stan and needed a confirmation:  
For a normal hierarchical case say  
Y ~ N(theta, sigma^2) with theta unknown and sigma^2 unknown.  
We know that the sample mean or Y-bar ~ N(theta, sigma^2 /n) where n is the sample size and

(n-1)\*S^2 / sigma^2 ~ Chi-square (degrees of freedom=(n-1)). With S^2 being the sample variance.

With a variable transformation we know under the normal case, S^2 ~ gamma((n-1)/2 , 2\*sigma^2/(n-1))  
How do you specify this in your model statement using “target +=” .  
Here was my implementation:

model  
{

target += normal\_lpdf(y\_bar|theta, sqrt(sigma2_inv(n)));  
target += gamma\_lpdf(S|n_.5-.5,(n-1)_inv(2.0_sigma2));

}

Further wanted to ask if I used the inv function correctly here or if there could be an additional function for specifying normal sufficient statistics(I saw a deprecated version). Last but not least, does it matter if you use an array or vector for later calculations of a parameter?  
Any help would be appreciated. For full stan code see below:

hb ← ’  
data  
{  
//Number of strata  
int\<lower=2\> K;

vector\<lower=0\>[K] n\_e; //sample sizes  
vector[K] y\_bar\_e; //sample means  
vector\<lower=0\> [K]S\_e; // sample variances

//pre-specified constants  
real k\_0;  
real b\_01;  
real b\_02;  
real b\_03;  
real b\_04;  
real\<lower=0\> v\_0;  
}

parameters  
{  
//level 1  
vector[K] theta\_ek;  
real\<lower=0.001\> sigma2\_e;

//level 2  
vector[K] mu\_k;  
vector[K] tau2\_k\_i;

//level 3  
real mu;  
real\<lower=0\> phi2\_i;  
real\<lower=0\> eta2\_i;  
}

model  
{  
// priors  
eta2\_i ~ gamma(b\_03, b\_04);  
phi2\_i ~ gamma(b\_01, b\_02);  
mu\_k ~ normal(2,sqrt(inv(phi2\_i)));  
tau2\_k\_i ~ gamma(v\_0_0.5,v\_0_0.5\*inv(eta2\_i));

theta\_ek ~ normal(mu\_k,sqrt(inv(tau2\_k\_i)));  
sigma2\_e ~ inv\_gamma(.01,.01);

// likelihood  
for (i in 1:K)  
{  
// Sample Means  
target += normal\_lpdf(y\_bar\_e[i]|theta\_ek[i], sqrt(sigma2\_e_inv(n\_e[i])));  
//Sample Variances  
target += gamma\_lpdf(S\_e[i]|n\_e[i]_.5-.5,(n\_e[i]-1)\* inv(2.0\*sigma2\_e));

}  
}  
’

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [June 24, 2021, 1:09pm UTC](https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171/2 "2021-06-24T13:09:30Z")

</div>

See [this post](https://discourse.mc-stan.org/t/psa-using-stantargets-for-sbc-esque-checks-big-speedups-for-big-gaussian-likelihoods/21425). And maybe [this one](https://discourse.mc-stan.org/t/looking-for-help-on-a-few-jacobians/22779) too.

---

<div class="post-metadata">

**Author:** ![eric11230](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric11230/32/11967_2.png) [@eric11230](https://discourse.mc-stan.org/u/eric11230)\
**Post date:** [June 24, 2021, 4:47pm UTC](https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171/3 "2021-06-24T16:47:31Z")

</div>

Thank you Mike! This helped. Do you know why I can not use target += as coded now?  
The model runs. Could also just not be a good model, but convergence is the issue and I think it’s because of the target +=.

---

<div class="post-metadata">

**Author:** ![eric11230](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric11230/32/11967_2.png) [@eric11230](https://discourse.mc-stan.org/u/eric11230)\
**Post date:** [June 8, 2022, 8:01pm UTC](https://discourse.mc-stan.org/t/sufficient-statistics-for-normal-case/23171/4 "2022-06-08T20:01:15Z")

</div>

I’ve re-explored this and can confirm that you can use sufficient statistics for the target += command.  
Below script 1 uses the normal likelihood and script 2 uses the normal sufficient statistics. Hope this can clear any doubts and help others better understand the target += command. By running both examples you will see the results agree.  
Example taken from: [Link](https://noah-padgett.github.io/Bayesian-Psychometric-Modeling/stan-model-for-mean-and-variance-unknown.html).

# Script 1

model\_normal ← ’  
data {  
int N;  
real x[N];  
real mu0;  
real sigma0;  
real alpha0;  
real beta0;  
}

parameters {  
real mu;  
real\<lower=0\> sigma\_2;  
}

model {  
x ~ normal(mu, sqrt(sigma\_2));  
mu ~ normal(mu0, sigma0);  
sigma\_2 ~ inv\_gamma(alpha0, beta0);  
}

’  
mydata ← list(  
N = 10,  
x=c(91, 85, 72, 87, 71, 77, 88, 94, 84, 92),  
mu0 = 75,  
sigma0 = 50,  
alpha0 = 5,  
beta0 = 150  
)

fit ← stan(  
model\_code = model\_normal,  
data = mydata,  
warmup = 1000,  
iter = 5000  
)

# Script 2

model\_normal1 ← ’  
data {  
int N;  
real xbar;  
real s2;  
real mu0;  
real sigma0;  
real alpha0;  
real beta0;  
}

parameters {  
real mu;  
real\<lower=0\> sigma\_2;  
}

model {

target += normal\_lpdf(xbar|mu, sqrt(sigma\_2/(N_1.0)));  
target += gamma\_lpdf(s2|N_.5-.5, N\*.5_inv(sigma\_2)-inv(2_sigma\_2));

mu ~ normal(mu0, sigma0);  
sigma\_2 ~ inv\_gamma(alpha0, beta0);  
}

’

real xbar;  
real s2;  
mydata ← list(  
N = 10,  
xbar=mean(c(91, 85, 72, 87, 71, 77, 88, 94, 84, 92)),  
s2=var(c(91, 85, 72, 87, 71, 77, 88, 94, 84, 92)),  
mu0 = 75,  
sigma0 = 50,  
alpha0 = 5,  
beta0 = 150  
)

fit ← stan(  
model\_code = model\_normal1,  
data = mydata,  
warmup = 1000,  
iter = 5000  
)