dev/monopoly_test.R

In ffeng23/ADASPR: An R package to characterize the distributed affinities in a multi-ligand system

library(MonoPoly)

#calculate the polynomial
beta <- c(1,2,1)
x <- 0:10
evalPol(x, beta)
str(evalPol(x, beta))
x <- cbind(0:10, 10:0)
evalPol(x, beta)
str(evalPol(x, beta))
x <- data.frame(x=0:10, y=10:0)
evalPol(x, beta)
str(evalPol(x, beta))

#fit the polynomial
fit <- monpol(y~x, w0)
ismonotone(fit)
beta <- c(1,0,2) ## the polynomial 1 + 2*x^2
ismonotone(beta)
ismonotone(beta, a=0)
ismonotone(beta, b=0)

###-------------------------------exp(x)
#testing the monotone fitting on 1-exp(x)
##R code to empirical verify exponential expansion
##
## y=exp(kx) <=> y=1+kx+(kx)^2/2!+(kx)^3/3!+(kx)^4/4!+(kx)^5/5!+....
##

x<-c(0:100)
degree<-120
fc<-factorial(c(1:degree))
k<- 1
xpoly<-poly(x*k,degree=degree,raw=T)

#####=======>note about poly
##poly is a function to calculate the polynomial format for the input
###in this case we have x =c(1,2,3,4,5,6,7,8,9,10). you can think of this 
##as 10 different concentration
# by calling poly with raw=TRUE, we got 
##[1] 1^1, 1^2,1^3, 1^4   #for the first conc=1 with degree up to 4
#@[2] 2^1, 2^2, 2^3, 2^4   #for the scond conc=2 with degree up to 4
##[3] 3^1,3^2, 3^3, 3^4  #for the third conc=3 with degree up to 4
##......
##......
#to get the full polynomial expression, we need to add x^0 terms to the matrix
#by calling cbind(xpoly, rep(1,dim(xpoly)[1]))


y<-exp(k*x)

y_poly<-rep(0,length(x))
for(i in c(1:length(x)))
{
y_poly[i]<-sum(xpoly[i,]/fc)+1
}


plot(x,y, log="", main="fitting exponential with polynomial")

lines(x,y_poly, col=2, lty=2 )
legend(5,0.8, c("exponential exp(-x)","polynomial with degree of 200"), col=c(1,2), pch=c(1,-1), lty=c(-1,2))

#now try to fit the polynomail
#x<-k*x 
mp<-monpol(y~x,degree=99,weights=1/exp(x))
xpoly<-poly(x,degree=99,raw=T)
xpoly<-cbind(x^0,xpoly)

#check for the 
b<-mp$beta;
yMp<-apply(xpoly,1,function(x){sum(x*b)})
lines(x,yMp, col=2,lty=2);