Hello all,
I am trying to use WinBUGS to estimate the parameters of a nonlinear
function, from 'data' with known Gaussian errors. (These 'data' are
simulated). The function is of the form:
f(x) = aa + bb*sin(cc*x^2 + dd*x + ee)
aa = 2.0
bb = 3.0
cc = 1.5
dd = -4.0
ee = 3.0
My 'data' consists of 30 (x,y) pairs. To generate each pair, I use the
following procedure:
1. choose x[i] from a uniform distribution on the range (0,5)
2. set y[i] = f(x[i]) + e[i]
where e[i] is selected from a Gaussian distribution with mean 0 and
standard deviation 0.2
I then use WinBUGS, through R, to estimate the parameters of the
function f(x) from the data.
My BUGS model is expressed in the file model_1.txt:
----------------%<--------------------------------------------------
# File model_1.txt
#
var x1[NDATA], y1[NDATA], yfit[NDATA], aa, bb, cc, dd, ee, tau, sigma
model
{
for (i in 1:NDATA)
{
yfit[i] <- aa + bb * sin(cc*x1[i]*x1[i] +dd*x1[i] + ee);
y1[i] ~ dnorm(yfit[i], tau);
}
tau <- 25.0
aa ~ dnorm(2.0, 0.01);
bb ~ dnorm(3.0, 0.01);
cc ~ dnorm(1.5, 0.01);
dd ~ dnorm(-4.0, 0.01);
ee ~ dnorm(3.0, 0.01);
}
----------------%<--------------------------------------------------
I run the calculation from R, using the R script model_1.R, below.
The variables NDATA, x1 and y1 contain the number of points (30), and
the x and y values of my simulated data.
----------------%<--------------------------------------------------
# File model_1.
#
data <- list("NDATA", "x1", "y1")
inits <- function()
{
aa = rnorm(1, 2.0, 5.0)
bb = rnorm(1, 3.0, 5.0)
cc = rnorm(1, 1.5, 5.0)
dd = rnorm(1, -4.0, 5.0)
ee = rnorm(1, 3.0, 5.0)
}
parameters = c("aa", "bb", "cc", "dd", "ee")
model.file = "model_1.txt"
fit1.sim <- bugs( data
, inits=NULL
, parameters
, model.file
, n.chains=5
, n.iter=100*1000
, coda=TRUE
, debug=TRUE
, DIC=TRUE
)
cda1 <- read.bugs(fit1.sim)
----------------%<--------------------------------------------------
Even with what I think is a very large number of iterations (100,000 --
which takes about 30 minutes on my laptop), it seems that the chains
have not converged. The resulting posterior distributions for each of
the five parameters is multimodal, and the trace plots (produced by
plot(cda1) in R) show no hint of convergence.
I do not understand what I might be doing wrong in my use of WinBUGS,
and would appreciate any assistance someone with more expertise might
provide.
best regards,
Dr. Marc Paterno
Fermi National Accelerator Laboratory
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|