Friday, October 17, 2014

NLS model in R

names(data) <- c("freq", "imp", "uuser", "click", "action");
data$CTR <- data$click/data$imp;
data$CVR <- data$action/data$imp;

data <- data[1:41,];
beta0 <- sum(data$action)/sum(data$imp);
beta1 <- max(data$action/data$imp);
beta2 <- -5/log(1/2);

decay <- nls(formula = CVR ~ b0 + b1 *exp(-freq/b2),
                 data=data,
                 start=c(b0=beta0, b1=beta1, b2=beta2),
                 trace = T)
summary(decay)
##Residual standard error: 5.08e-06 on 137648 degrees of freedom
## the coefficients only:
coef(decay)
## including their SE, etc:
coef(summary(decay))
plot(data$freq,data$CVR, main = "nls(*), data, true function and fit, n=100")
#curve(beta0 + beta1 * exp(-x / beta2), col = 4, add = TRUE)
lines(data$freq, predict(decay), col = 2)

result <- data.frame(data$freq,data$CVR,predict(decay));

ggplot(data, aes(y=data$CVR, x=data$freq)) +
  geom_point(stat = "identity", fill = "red", size = 3) +
#  stat_smooth() +
  geom_line(aes(y=predict(decay), x=data$freq), color = 'red', size = 1, lty = 2) +
  ggtitle("CVR by Frequency") +
  ylab("CVR") + xlab("Frequency in one day") +
  theme(plot.title = element_text(face = "bold", size = 20)) +
  theme(axis.text.x = element_text(face = "bold", size = 16)) +
  theme(axis.text.y = element_text(face = "bold", size = 16)) +
  theme(axis.title.x = element_text(face = "bold", size = 16)) +
  theme(axis.title.y = element_text(face = "bold", size = 16, angle = 90)) +
  theme(legend.position = "top") +
  theme(legend.key = element_rect(colour = NA)) +
  theme(legend.title = element_blank()) +
  scale_x_discrete(breaks = as.character(seq(0, 30, by = 2)))

##starting value
negexp.sv <- function(x,y){
  mx <- mean(x)
  x1<-x-mx
  x2<-((x-mx)^2)/2
  b <- as.vector(lm(y~x1+x2)$coef)
  b2 <- b[2]/b[3]
  b<-as.vector(lm(y~exp(-x/b2))$coef)
  parms <- cbind(b[1],b[2],b2)
}
b.start <- negexp.sv(data$freq, data$CVR)
##specify gradient
expn <- function(b0,b1,b2,x){
  temp <- exp(-x/b2)
  model.func <- b0+b1*temp
  D <- cbind(1, temp, (b1*x*temp)/b2^2)
  dimnames(D) <- list(NULL, c("b0","b1","b2"))
  attr(model.func, "gradient") <- D
  model.func
}

decay <- nls(formula = CVR ~ expn(b0, b1, b2, frequency),
             data=data,
##             start=c(b0=b.start[1], b1=b.start[2], b2=b.start[3]),
             start=c(b0=beta0, b1=beta1, b2=beta2),
             trace = T)
summary(decay)
##Residual standard error: 5.08e-06 on 137648 degrees of freedom


for (i in seq(300,1000,100)){
  print(i);
  data1 <- data[1:i,];
  decay <- nls(formula = CVR ~ expn(b0, b1, b2, frequency),
               data=data1,
##               start=c(b0=b.start[1], b1=b.start[2], b2=b.start[3]),
               start=c(b0=beta0, b1=beta1, b2=beta2),
               trace = F)
  print(coef(summary(decay)))


}

No comments:

Post a Comment