Feature request: Linear combinations of hyperparameters

[GoogleCodeExporter]

An option to compute the marginals of linear combinations of hyperparameters is 
necessary to correctly work with the variance of the SPDE models.  With this 
feature, we can modify the 'plot' and 'summary' functions to display the 
posteriors in terms of the intuitive parameterisation.

Original issue reported on code.google.com by dp.simp...@gmail.com on 10 May 2011 at 2:19

[GoogleCodeExporter]

Do you mean linear combinations of only hyperparameters, or  joint linear 
combinations of the hyperparameters and the latent field ? with hyperparameters 
you mean hyperparameters in the internal or external scale ? (like 
log(precision) and precision). 

Original comment by havard.rue on 16 May 2011 at 4:36

[GoogleCodeExporter]

The main point is to get the marginal for the variance of the SPDE model, which 
requires a linear combination of the the log(kappa) and log(tau) marginals.  It 
would also be good to be able to get posteriors for the sums of basis functions 
in the nonstationary case (this would need to be at a massive number of points, 
hence it needs to work for a large number of linear combinations).

Original comment by dp.simp...@gmail.com on 16 May 2011 at 6:21

[GoogleCodeExporter]

UserFuncion0/1 currently gives the marginals for non-stationary log(tau) and 
log(kappa).  Whet is needed is a generalised version where the user specifies 
the "basis functions"="linear combinations of the parameters" separately from 
the ones used for tau and kappa.  This will also later tie in with the planned 
reparameterisation of tau and kappa into a more generalised version.

Original comment by finn.lin...@gmail.com on 17 May 2011 at 7:28

[GoogleCodeExporter]

Clarification:  The difference to the current UserFunction is that the theta.T 
and theta.K parameters need to enter into the same linear combinations.  Until 
the new model parameterisation is implemented, we don't need a fully general 
solution, but can get away with just extending the "userfunctions" with one 
more lincomb-field (the range marginals can be obtained by a univariate 
nonlinear transformation of the log(kappa^2) field): The field variances are 
obtained from a univariate nonlinear transformation of 2*log(tau)+log(kappa^2)

Original comment by finn.lin...@gmail.com on 17 May 2011 at 7:35

[GoogleCodeExporter]

I also wanted this and after some discussion Sara Martino wrote a function to 
sample from the joint distribution of theta, e.g.,

fit <- inla(....)
theta <- theta.sample(fit)

which gives a sample of the hyerparameters in the internal scale (e.g. log 
precisions) and from those one can compute linear combinations like we would 
have MCMC samples. I do not know how accurate that is.

Original comment by gregor.g...@gmail.com on 3 Jun 2011 at 10:23

[GoogleCodeExporter]

For hyperparameters, I tend to think that a such solution would be sufficient. 
First because it more common to look at non-linear functions of the 
hyperparameters, like the sum of two variances, which would look like

    exp(-theta1) + exp(-theat2)

using the internal representation. the dimension is also small, so a sampling 
based approach which cover all cases would be ok until a better approach is 
available. The only issue is that the mapping between the real scale (f.ex 
precision) and the internal scale (f.ex log precision) is not there in the 
output, so the user must know what (s)he is doing. it might be possible to pass 
these to.theta() from.theta() functions in models.R through Model.ini and then 
back again to the results, as plain text, which can then be evaluated when 
reading the results. hence, we could provide samples from theta in the user 
scale as well, automatically. 

H

-- 
H�vard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533    URL  : http://www.math.ntnu.no/~hrue  
Fax  : +47-7359-3524    Email: havard.rue@math.ntnu.no

This message was created in a Microsoft-free computing environment.

Original comment by havard.rue on 4 Jun 2011 at 7:43

[GoogleCodeExporter]

The user will (and should be allowed to) assume that to.theta and from.theta 
are valid for the output as well as the input, for the cases where they are 
well-defined, so no need to pass them back.

It would be useful to have wrapper functions however, to automatically 
translate all thetas into their documented "untransformed" states.  This will 
in particular be true for the coming general spde model parameterisations, 
where the internal representation will be variable, depending on options to 
inla.spde.create.

Helper functions for combining theta-samples with non-linear (and 
multiple-valued) functions with the transformations will then be needed to more 
easily obtain quantiles etc.

Original comment by finn.lin...@gmail.com on 4 Jun 2011 at 8:03

[GoogleCodeExporter]

I just added (and rebuildt the package) the function

inla.sample.hyper(n, result)

which creates samples from theta in user or internal scale using the CCD 
interpolant. (i.e. the Gaussian with adjusted stdevs in each direction).   

The function will fail if used with 'group' if samples from theta in user-scale 
is requested (default), as these to.theta from.theta functions require number 
of groups as well as argument. 

this is a test-implementation and its currently undocumented. 

let me know what you think or if you have suggestions to improvements.

H

Original comment by havard.rue on 8 Jun 2011 at 8:52

[GoogleCodeExporter]

if(intern) then use "Precision" instead of "Log precision" in colum names

Original comment by gregor.g...@gmail.com on 21 Nov 2011 at 12:02

[GoogleCodeExporter]

Thanks for this request. As the function does not know the 'names' of the 
hyperparameters, it adds ``in user-scale'' to all the names when intern is 
TRUE. Doing it correct, its kind of messy... (I blame myself on this issue!)

I'll put it on the list. 

H

> r = inla(y~1+f(i),data = data.frame(y=1,i=1))
> inla.sample.hyper(n=4,r,intern=TRUE)
         Log precision for the Gaussian observations Log precision for i
sample-1                                 9.490503285         9.809773793
sample-2                                 8.924006093         8.384219807
sample-3                                10.645525978         9.477518832
sample-4                                 7.379825217         8.569078704
> inla.sample.hyper(n=4,r,intern=FALSE)
         Log precision for the Gaussian observations in user-scale
sample-1                                              31042.457556
sample-2                                              24334.405346
sample-3                                               7721.941180
sample-4                                               4320.985485
         Log precision for i in user-scale
sample-1                       10486.63470
sample-2                       97086.05944
sample-3                       12353.64269
sample-4                       25774.09746

Original comment by havard.rue on 21 Nov 2011 at 7:55

[GoogleCodeExporter]

... and I suggest to rename this function to inla.hyperpar.sample in accordance 
with inla.hyperpar.

gg

Original comment by gregor.g...@gmail.com on 21 Nov 2011 at 11:16

[GoogleCodeExporter]

It is of vital importance that the names be machine-interpretable, which in 
this case means that they must be constructed in a well-specified way from the 
canonical parameter names, e.g.
inla.models()$ar1$hyper$theta1$name
or
inla.models()$ar1$hyper$theta1$short.name

At the very least, there needs to be an option requesting such behaviour (e.g. 
choosing between "short", "long", and "human").  This goes for _all_ output 
from inla...

Original comment by finn.lin...@gmail.com on 21 Nov 2011 at 11:28

[GoogleCodeExporter]

Function is now renamed to inla.hyperpar.sample().

Original comment by havard.rue on 21 Nov 2011 at 7:42