## Data
y=iris$Petal.Width
x=iris$Petal.Length

## Set Basis Knots
knot=seq(2,6,1)

K=length(knot)
n1=rep(1,length(x))
n2=x
n3=d(knot[1],knot[K],x)-d(knot[(K-1)],knot[K],x)
n4=d(knot[2],knot[K],x)-d(knot[(K-1)],knot[K],x)
n5=d(knot[3],knot[K],x)-d(knot[(K-1)],knot[K],x)

## d_K Function
d=function(ck,cK,x){
  
  dkx={0}
   for ( i in 1:length(x) ){
     
     fir=ifelse(x[i]>ck,(x[i]-ck)^3,0)
     sec=ifelse(x[i]>cK,(x[i]-cK)^3,0)
     dkx[i]=(fir-sec)/(cK-ck)
     
   }  
  return(dkx)
  
}

## Basis
N=cbind(n1,n2,n3,n4,n5)

## Get beta_hat from OLS
bhat = solve(t(N)%*%N)%*%t(N)%*%y

## Give any x in domain for x 
prex=seq(1.5,6.5,0.5)

## Sub prex to the same basis
pren1=rep(1,length(prex))
pren2=prex
pren3=d(knot[1],knot[K],prex)-d(knot[(K-1)],knot[K],prex)
pren4=d(knot[2],knot[K],prex)-d(knot[(K-1)],knot[K],prex)
pren5=d(knot[3],knot[K],prex)-d(knot[(K-1)],knot[K],prex)
preN=cbind(pren1,pren2,pren3,pren4,pren5)

## pre = sum_i beta_hat * h(x_i)
pre=as.vector(t(bhat) %*% t(preN))

plot(x,y)
lines(x=prex,y=pre)