# Estimate The Rate At Which Standard Deviation Scales With An Independent Variable

I have an experiment in which I am taking measurements of a normally distributed variable $Y$,

$$Y \sim N(\mu,\sigma)$$

However, previous experiments have provided some evidence that the standard deviation $\sigma$ is an affine function of an independent variable $X$, i.e.

$$\sigma = a|X| + b$$

$$Y \sim N(\mu,a|X| + b)$$

I would like to estimate the parameters $a$ and $b$ by sampling $Y$ at multiple values of $X$. Additionally, due to experiment limitations I can only take a limited (roughly 30-40) number of samples of $Y$, and would prefer to sample at several values of $X$ for unrelated experimental reasons. Given these constraints, what methods are available to estimate $a$ and $b$?

Experiment Description

This is extra information, if you're interested in why I'm asking the above question. My experiment measures auditory and visual spatial perception. I have an experiment setup in which I can present either auditory or visual targets from different locations, $X$, and subjects indicate the perceived location of the target, $Y$. Both vision* and audition get less precise with increasing eccentricity (i.e. increasing $|X|$), which I model as $\sigma$ above. Ultimately, I'd like to estimate $a$ and $b$ for both vision and audition, so I know the precision of each sense across a range of locations in space. These estimates will be used to predict the relative weighting of visual and auditory targets when presented concurrently (similar to the theory of multisensory integration presented here: http://www.ncbi.nlm.nih.gov/pubmed/12868643).

*I know that this model is inaccurate for vision when comparing foveal to extrafoveal space, but my measurements are constrained solely to extrafoveal space, where this is a decent approximation.

• Interesting problem. It's likely the best solutions will take into account the reasons you are doing this experiment. What are your ultimate objectives? Prediction? Estimation of $\mu$, $a$, and/or $\sigma$? The more you can tell us about the purpose, the better the answers can be. – whuber Jan 9 '15 at 21:32
• Since the SD cannot be negative, it is unlikely to be a linear function of X. Your suggestion, a|X|, necessitates a narrower or wider V shape w/ a minimum at X=0, which seems a rather unnatural possibility to me. Are you sure this is right? – gung - Reinstate Monica Jan 9 '15 at 21:44
• Good point @gung, I'd inappropriately oversimplified my problem. It'd be more realistic to say that $\sigma$ is an affine function of $|X|$. I'll edit my question. – Adam Bosen Jan 9 '15 at 21:52
• @whuber The reason for wanting this is a bit involved, but I'll think about how to explain the experiment and add some more detail to my question soon. – Adam Bosen Jan 9 '15 at 21:56
• Do you have good reason, a-priori, to believe that X=0 represents the minimum SD, & that f(|X|) is monotonic? – gung - Reinstate Monica Jan 9 '15 at 22:03

In a case like yours, where you have a relatively simple, but "non-standard" generative model that you'd like to estimate parameters for, my first thought would be to use a Bayesian inference program like Stan. The description you've given would translate very cleanly to a Stan model.

Some example R code, using RStan (the R interface to Stan).

library(rstan)

model_code <- "
data {
int<lower=0> n; // number of observations
real y[n];
real x[n];
}
parameters {
real mu; // I've assumed mu is to be fit.
// Move this to the data section if you know the value of mu.
real<lower=0> a;
real<lower=0> b;
}
transformed parameters {
real sigma[n];
for (i in 1:n) {
sigma[i] <- a + b * fabs(x[i]);
}
}
model {
y ~ normal(mu, sigma);
}
"

# Let's generate some test data with known parameters

mu <- 0
a <- 2
b <- 1

n <- 30
x <- runif(n, -3, 3)
sigma <- a + b * abs(x)
y <- rnorm(n, mu, sigma)

# And now let's fit our model to those "observations"

fit <- stan(model_code=model_code,
data=list(n=n, x=x, y=y))

print(fit, pars=c("a", "b", "mu"), digits=1)


You'll get output that looks something like this (although your random numbers will probably be different to mine):

Inference for Stan model: model_code.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean  sd 2.5%  25% 50% 75% 97.5% n_eff Rhat
a   2.3       0 0.7  1.2  1.8 2.2 2.8   3.9  1091    1
b   0.9       0 0.5  0.1  0.6 0.9 1.2   1.9  1194    1
mu  0.1       0 0.6 -1.1 -0.3 0.1 0.5   1.4  1262    1

Samples were drawn using NUTS(diag_e) at Thu Jan 22 14:26:16 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).


The model has converged well (Rhat=1), and the effective sample size (n_eff) is reasonably large in all cases, so on a technical level the model is well-behaved. The best estimates of $a$, $b$ and $\mu$ (in the mean column) are also fairly close to what was provided.

• Oh, I like this! I hadn't heard of Stan before, thanks for the reference. I was initially hoping for an analytic solution, but given the lack of responses I doubt one exists. I'm inclined to believe your answer is the best approach to this problem. – Adam Bosen Jan 22 '15 at 19:49
• It wouldn't completely shock me if an analytic solution existed, but I'd certainly be a little bit surprised. The strength of using something like Stan is that it's very easy to make changes to your model - an analytic solution would probably be very strongly constrained to the model as given. – Martin O'Leary Jan 23 '15 at 12:09

You cannot expect closed formulas, but you can still write down the likelihood function and maximize it numerically. Your model is $$\newcommand{\dist}{\sim} Y \dist N(\mu, a|x|+b)$$ Then the loglikelihoodfunction (apart from a term not depending on parameters) becomes $$l(\mu, a, b) = -\sum \ln(a|x_i|+b) -\frac12\sum\left(\frac{y_i-\mu}{a|x_i|+b}\right)^2$$ and that is easy to program and give to a numerical optimizer.

In R, we can do

make_lik  <-  function(x,y){
x  <-  abs(x)
function(par) {
mu <- par[1];a  <-  par[2];  b <-  par[3]
axpb <-  a*x+b
-sum(log(axpb)) -0.5*sum( ((y-mu)/axpb)^2 )
}
}


Then simulate some data:

> x <-  rep(c(2,4,6,8),10)
> x
[1] 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4
[39] 6 8
> a <- 1
> b<-  3
> sigma <-  a*x+b
> mu  <-  10
> y  <-  rnorm(40,mu, sd=sigma)


Then make the loglikelihood function:

> lik <-  make_lik(x,y)
> lik(c(10,1,3))
[1] -99.53438


Then optimize it:

> optim(c(9.5,1.2,3.1),fn=function(par)-lik(par))
$par [1] 9.275943 1.043019 2.392660$value
[1] 99.12962

$counts function gradient 136 NA$convergence
[1] 0

\$message
NULL