Finding a trend in variance

Question

I have written a little function that samples from a normal distribution with linearly changing variance

rnormltv <- function(n,mu,rsd) {
    vc <- vector(mode="numeric", length=n)
    i <- 1
    for (x in seq(from=rsd[1], to=rsd[2], length=n)) {
        vc[i]=rnorm(1,mu,x)
        i <- i+1
    }
    return(vc)
}

Suppose we obtain 1000 samples, with a decreasing variance from 9 to 1.

set.seed(1)
x <- rnormltv(1000,0,c(9,1))
plot(x)

The goal is to estimate the parameters 9 and 1 given this data, the mean of zero, assuming a normal distribution and a linear trend in variance. Conceptually I can see how this can be done with maximum likelihood, but how can it by done in practice, preferably in R ?

Note: This is not homework.

Answer 1

Here is a similar solution in R.

Please note that the code you provided generates random normal variables with a linear trend in sthe standard deviation, not in the variance. So I used this, you can easily modify the code for other kind of trends.

The only caveat is to compute directly the log likelihood, which will behave better with big sample, where the likelihood is very near to 0

set.seed(1)
x <- rnormltv(1000,0,c(9,1))
f <- function(a,b,x) sum( dnorm(x, mean=0, sd = a+(b-a)*seq(0,1,length=length(x)), log=TRUE) )

nlm( function(p) -f(p[1], p[2], x), p = c(8,1) )
$minimum
[1] 2923.871

$estimate
[1] 9.120767 1.126542

$gradient
[1] 4.048507e-05 2.179800e-05

$code
[1] 1

$iterations
[1] 9

So the estimate is $a = 9.12$ and $b=1.13$.

Please follow whuber advice to check the robustness of this estimation to your assumptions (zero mean, linear trend in standard deviation). If you need this for a particular kind of data analysis, always have a look on what other people did with the same kind of data.

Answer 2

You mentioned initially trying Maple. I did, too, and here's what I got (disclosure: I work for them):

# To generate a test sample:
decayingSample := proc(N :: nonnegint, a :: positive, b :: positive, $)
local i, sd, sample, result;
  result := Vector(N, 'datatype=float');
  sample := Statistics:-Sample(Normal(0, sd));
  for i to N do
    sd := a + (i-1)*(b-a)/(N-1);
    result[i] := sample(1)[1];
  end do;
  return result;
end proc:

# aa, bb are the values for our test sample; a and b will represent
# the unknown symbolic parameters to be determined.
NN, aa, bb := 1000, 9, 1:
ds := decayingSample(NN, aa, bb):
with(Statistics):
PointPlot(ds);

# Form the symbolic log likelihood for the ith entry.
ll := LogLikelihood(Normal(0, normal(a + (i-1)*(b-a)/(NN-1))),
                    x, samplesize=1) assuming 'real';
# Now add these guys for i = 1 .. NN:
expr := evalf(add('eval(ll, ['i' = i, x[1] = ds[i]])', i = 1..NN)):
# Create an initial guess for a and b:
initialGuess := [a = StandardDeviation(ds[1 .. iquo(NN, 2)]),
                 b = StandardDeviation(ds[iquo(NN, 2) .. NN])];
# Optionally:
# infolevel[Optimization] := 5: # for more detailed info
Optimization:-Maximize(expr, {a >= 0, b >= 0}, 'initialpoint' = initialGuess);

This lead to an estimate of $a = 8.72, b = 1.07$; trying it on a second sample gave $a = 9.20, b = 0.94$.