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.
 A: 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.
A: 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

