[R] Constrained cubic smoothing spline

2013-03-05 Thread Victor hyk
Hello everone,
   Anyone who knows how to force a cubic smoothing spline to pass 
through a particular point?
   I found on website  someone said that we can use "cobs package" to 
force the spline pass through certain points or impose shape   
constraints (increasing, decreasing). However,  this package is using  B-spline 
and can only do linear and quadratic B-spline.
   In my research, I need to force a cubic smoothing spline to pass a 
point. 
   Thanks!

  Victor

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] About basic logical operators

2013-03-05 Thread Victor hyk
Hello everyone,
  I have a basic question regarding logical operators.
> x<-seq(-1,1,by=0.02)
> x
  [1] -1.00 -0.98 -0.96 -0.94 -0.92 -0.90 -0.88 -0.86 -0.84 -0.82 -0.80 -0.78
 [13] -0.76 -0.74 -0.72 -0.70 -0.68 -0.66 -0.64 -0.62 -0.60 -0.58 -0.56 -0.54
 [25] -0.52 -0.50 -0.48 -0.46 -0.44 -0.42 -0.40 -0.38 -0.36 -0.34 -0.32 -0.30
 [37] -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06
 [49] -0.04 -0.02  0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18
 [61]  0.20  0.22  0.24  0.26  0.28  0.30  0.32  0.34  0.36  0.38  0.40  0.42
 [73]  0.44  0.46  0.48  0.50  0.52  0.54  0.56  0.58  0.60  0.62  0.64  0.66
 [85]  0.68  0.70  0.72  0.74  0.76  0.78  0.80  0.82  0.84  0.86  0.88  0.90
 [97]  0.92  0.94  0.96  0.98  1.00
> x[x<=0.02]
 [1] -1.00 -0.98 -0.96 -0.94 -0.92 -0.90 -0.88 -0.86 -0.84 -0.82 -0.80 -0.78
[13] -0.76 -0.74 -0.72 -0.70 -0.68 -0.66 -0.64 -0.62 -0.60 -0.58 -0.56 -0.54
[25] -0.52 -0.50 -0.48 -0.46 -0.44 -0.42 -0.40 -0.38 -0.36 -0.34 -0.32 -0.30
[37] -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06
[49] -0.04 -0.02  0.00
> x[x<0.2]
 [1] -1.00 -0.98 -0.96 -0.94 -0.92 -0.90 -0.88 -0.86 -0.84 -0.82 -0.80 -0.78
[13] -0.76 -0.74 -0.72 -0.70 -0.68 -0.66 -0.64 -0.62 -0.60 -0.58 -0.56 -0.54
[25] -0.52 -0.50 -0.48 -0.46 -0.44 -0.42 -0.40 -0.38 -0.36 -0.34 -0.32 -0.30
[37] -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06
[49] -0.04 -0.02  0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18
[61]  0.20
> 
 Why does "x[x<=0.02]" return  no 0.02 but "x[x<0.2]" return a subsample 
with 0.02?
 Anyone who can tell me why?
 Thanks!

 Victor

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Use pcls in "mgcv" package to achieve constrained cubic spline

2013-03-11 Thread Victor hyk
Hello everyone,
 Dr. wood told me that I can adapting his example to force cubic spline 
to pass through certain point.
 I still have no idea how to achieve this. Suppose we want to force the 
cubic spline to pass (1,1), how can 
I achieve this by adapting the following code?

# Penalized example: monotonic penalized regression spline .
# Generate data from a monotonic truth.
x<-runif(100)*4-1;x<-sort(x);
f<-exp(4*x)/(1+exp(4*x));y<-f+rnorm(100)*0.1;plot(x,y)
dat<-data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug<-gam(y~s(x,k=10,bs="cr"));lines(x,fitted(f.ug))
# Create Design matrix, constraints etc. for monotonic
 spline
sm<-smoothCon(s(x,k=10,bs="cr"),dat,knots=NULL)[[1]]
F<-mono.con(sm$xp); # get constraints
G<-list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=sm$xp,y=y,w=y*0+1)
G$Ain<-F$A;G$bin<-F$b;G$S<-sm$S;G$off<-0
p<-pcls(G); # fit spline (using s.p. from unconstrained fit)
fv<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv,col=2)
lines(x,f,col="blue") 
   Thanks a lot!!

    Victor
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Impose constraint on first order derivative at a point for cubic smoothing spline

2013-11-01 Thread Victor hyk
Hello, 
      Dr. Simon Wood told me how to force a cubic spline passing through a
point. The code is as following. Anyone  who knows how I can change the code
to force the first derivative to be certain value. For example, the first
derivative of the constrained cubic spline equals 2 at point (0, 0.6). 
      I really appreciate your help! 
      Thanks! 
        
      Best      
      Victor 
      
      Here is the initial reply and code provided by Dr. Simon Wood: 

"Actually, you might as well use "gam" directly for this (has the 
advantage that the smoothing parameter will be chosen correctly subject 
to the constraint). Here is some code. Key idea is to set the basis and 
penalty for the spline up first, apply the constraint, and then use gam 
to fit it..." 

best, 
Simon 

## Example constraining spline to pass through a 
## particular point (0,.6)... 

## Fake some data... 

library(mgcv) 
set.seed(0) 
n <- 100 
x <- runif(n)*4-1;x <- sort(x); 
f <- exp(4*x)/(1+exp(4*x));y <- f+rnorm(100)*0.1;plot(x,y) 
dat <- data.frame(x=x,y=y) 

## Create a spline basis and penalty, making sure there is a knot 
## at the constraint point, (0 here, but could be anywhere) 
knots <- data.frame(x=seq(-1,3,length=9)) ## create knots 
## set up smoother... 
sm <- smoothCon(s(x,k=9,bs="cr"),dat,knots=knots)[[1]] 

## 3rd parameter is value of spline at knot location 0, 
## set it to 0 by dropping... 
X <- sm$X[,-3]        ## spline basis 
S <- sm$S[[1]][-3,-3] ## spline penalty 
off <- y*0 + .6      ## offset term to force curve through (0, .6) 

## fit spline constrained through (0, .6)... 
b <- gam(y ~ X - 1 + offset(off),paraPen=list(X=list(S))) 
lines(x,predict(b)) 

## compare to unconstrained fit... 

b.u <- gam(y ~ s(x,k=9),data=dat,knots=knots) 
lines(x,predict(b.u))  

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.