Description Usage Arguments Details Value Note References See Also Examples
View source: R/curve.polynomial.rjmcmc.R
Fit a variety of curves by a sequence of piecewise polynomials. The number and location of knots are determined by the Reversible Jump MCMC method.
1 2 3 | curve.polynomial.rjmcmc(y, x, lambda, l, l0, c=0.4,
gamma.shape=1e-3, gamma.rate=1e-3,
maxit=10000, err=1e-8, verbose=TRUE)
|
y |
a vector of values of a response variable. |
x |
a vector of values of the corresponding explanatory variable. |
lambda |
the parameter of the Poisson prior for the number of knots. |
l |
the order of polynomials. |
l0 |
the degree of continuity at the knots. |
c |
the parameter controlling the rate of changing dimension. The default value is 0.4. |
gamma.shape |
the parameter shape of the gamma prior for the error precision. The default value is 1e-3. |
gamma.rate |
the parameter rate of the gamma prior for the error precision. The default value is 1e-3. |
maxit |
the maximum number of iterations. The default value is 10000. |
err |
the iteration stops when consecutive difference in percentage of mean-squared error reaches this bound. The default value is 1e-8. |
verbose |
logical. If |
The method is based on Denison et al. (1998). It can be used to fit
both smooth and unsmooth curves. When both l0
and l
are set to 3, it fits curves using cubic spline.
knots.save |
a list of location of knots. One component per iteration. |
betas.save |
a list of estimates of regression parameters βs. One component per iteration. |
fitted.save |
a matrix of fitted values. One column per iteration. |
sigma2.save |
a vector of estimate of error variance. One component per iteration. |
The factor \frac{1}{√{n}} was added in the likelihood ratio to penalize the ratio for dimensionality as suggested in Dimatteo et al. (2001).
Denison, D. G. T., Mallick, B. K., Smith, A. F. M.(1998) Automatic Bayesian Curve Fitting JRSSB vol. 60, no. 2 333-350
Dimatteo, I., Genovese, C. R., Kass, R. E.(2001) Bayesian Curve-fitting with Free-knot Splines Biometrika vol. 88, no. 4 1055-1071
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | ## Not run:
#Example 1: smooth curve
#example 3.1. (b) in Denison et al.(1998)
x <- runif(200)
xx <- -2 + 4*x
y.truth <- sin(2*xx) + 2*exp(-16*xx^2)
y <- y.truth + rnorm(200, mean=0, sd=0.3)
results <- curve.polynomial.rjmcmc(y, x, lambda=1, l=2, l0=1)
plot(sort(x), y.truth[order(x)], type="l")
lines(sort(x), rowMeans(results$fitted.save), col="red")
#Example 2: unsmooth curve
#blocks in Denison et al.(1998)
tj <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81)
hj <- c(4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2)
t <- seq(0, 1, len=2048)
Ktmtj <- outer(t, tj, function(a,b) ifelse(a-b > 0, 1, -1))
ft <- rowSums(Ktmtj %*% diag(hj))
x <- t
y <- ft + rnorm(2048, 0, 1)
results <- curve.polynomial.rjmcmc(y, x, lambda=5, l=2, l0=1)
plot(x, ft, type="s")
lines(x, rowMeans(results$fitted.save), col="red")
#Example 3: unsmooth curve
#bumps in Denison et al.(1998)
tj <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81)
hj <- c(4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)*10
wj <- c(0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005)
t <- seq(0, 1, len=2048)
ft <- rowSums((1 + abs(outer(t, tj ,"-") %*% diag(1/wj)))^(-4) %*% diag(hj))
y <- ft + rnorm(2048, 0, 1)
results <- curve.polynomial.rjmcmc(y, t, lambda=5, l=2, l0=1)
plot(t, ft, type="s")
lines(t, rowMeans(results$fitted.save), col="red")
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.