####Density Estimation


rec <- function(x) (abs(x) < 1) * 0.5
tri <- function(x) (abs(x) < 1) * (1 - abs(x))
gauss <- function(x) 1/sqrt(2*pi) * exp(-(x^2)/2)
x <- seq(from = -3, to = 3, by = 0.001)
plot(x, rec(x), type = "l", ylim = c(0,1), ylab = expression(K(x)),xlab="x")
lines(x, tri(x), lty = 2)
lines(x, gauss(x), lty = 3)
legend("topleft", legend = c("Rectangular", "Triangular", "Gaussian"), lty = 1:3, title = "kernel functions", bty = "n")

x <- c(0, 1, 1.1,1.2,1.3,1.4,1.5, 1.9, 2.8, 2.9, 3.5);x
n <- length(x)

bw.SJ(x)
plot(density(x,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="x")
lines(density(x, bw = "nrd"), col = 2)
lines(density(x, bw = "ucv"), col = 3)
lines(density(x, bw = "bcv"), col = 4)
lines(density(x, bw = "SJ-ste"), col = 5)
lines(density(x, bw = "SJ-dpi"), col = 6)
legend("topleft",
       legend = c("nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi"),
       col = 1:6, lty = 1)


####Fitting density function with package#####
a=file.choose()
CMB=read.table(a,header=TRUE)
names(CMB)
head(CMB)
plot(CMB)
multipole=CMB[,1]
power=CMB[,2]

bw.SJ(power)
plot(density(power,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="x")
bw.nrd0(multipole)
bw.SJ(multipole)
plot(density(multipole,bw="nrd0",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="x")
rug(multipole)
boxplot(multipole,ylab="Score")

plot(density(multipole, n = 1000),xlab="")
rug(multipole)





require(graphics)

plot(density(power))
rug(power)
lines(density(power, bw = "nrd"), col = 2)
lines(density(power, bw = "ucv"), col = 3)
lines(density(power, bw = "bcv"), col = 4)
lines(density(power, bw = "SJ-ste"), col = 5)
lines(density(power, bw = "SJ-dpi"), col = 6)
legend(55, 0.035,
       legend = c("nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi"),
       col = 1:6, lty = 1)

#######################################Non-parametric Regression
#################
####exmapel 1: ####
#################
a=file.choose()
CMB=read.table(a,header=TRUE)
names(CMB)
head(CMB)
plot(CMB)
multipole=CMB[,1]
power=CMB[,2]
plot(multipole,power,xlab="multipole",ylab = "power")
abline(v=400,col="red")
plot(multipole[1:400],power[1:400],xlab="multipole",ylab = "power")
#################
####exmapel 2####
#################
?ksmooth
layout(matrix(1:4,nrow=2,ncol = 2))
plot(multipole[1:400],power[1:400],main="h=1",pch=1,ylab="Power",xlab="Multipole")
lines(ksmooth(CMB$multipole,CMB$power,bandwidth = 1),lwd=2,col="red")
plot(multipole[1:400],power[1:400],main="h=10",pch=1,ylab="Power",xlab="Multipole")
lines(ksmooth(CMB$multipole,CMB$power,bandwidth = 10),lwd=2,col="red")
plot(multipole[1:400],power[1:400],main="h=50",pch=1,ylab="Power",xlab="Multipole")
lines(ksmooth(CMB$multipole,CMB$power,bandwidth = 50),lwd=2,col="red")
plot(multipole[1:400],power[1:400],main="h=200",pch=1,ylab="Power",xlab="Multipole")
lines(ksmooth(CMB$multipole,CMB$power,bandwidth =200),lwd=2,col="red")
#################
####exmapel 3####
#################
library(sm)
?hcv
hn=hcv(CMB$multipole,CMB$power, display="lines",ngrid=260)
hn
sm.regression(multipole,power, h=hn,display = "image", 
       theta = 120)
par(mfrow = c(2,1))

plot(multipole,power,ylab="power",xlab="multipole",xlim=c(0,900),ylim=c(0,7000),pch=".")
lines(ksmooth(multipole, power, "normal", bandwidth=hn),lwd=2,col="green")
#################
###exmapel 4####
#################
library(SemiPar)
library(locfit)
?lidar
data(lidar)
head(lidar)
names(lidar)

par(mfrow=c(2,2))
l=lidar[,1]
r=lidar[,2]
hn=hcv(l,r, display="lines",ngrid=221)
hn
lcv(l,r)
plot(l,r, main="local linear regression",col="pink",pch=20)
lines(locfit(lidar$logratio ~ lidar$range,alpha=c(0,37)),lwd=2,col="green")
library(KernSmooth)
lines(locpoly(l, r, bandwidth=15, degree=4), col="blue",lwd=2)

ei=resid(locfit(lidar$logratio ~ lidar$range))
plot(lidar$range,ei,ylim=c(-0.4,0.4),xlim=c(350,800),xlab="range",ylab="Residuals")


#################
####exmapel 5####
#################

library(locfit)
hn=hcv(CMB$multipole, CMB$power, display="lines",ngrid=50)
hn
layout(matrix(1:4,nrow=2,ncol = 2))

plot(multipole[1:400],power[1:400], main="Locally polynomials with p=0",ylab="",xlab="")

lines(locpoly(multipole[1:400],power[1:400],bandwidth =hn,deg=0),lwd=3,col="red")

plot(multipole[1:400],power[1:400], main="Locally polynomials with p=1",ylab="",xlab="")

lines(locpoly(multipole[1:400],power[1:400],bandwidth =hn,deg=1),lwd=3,col="red")

plot(multipole[1:40],power[1:40], main="Locally polynomials with p=0",ylab="",xlab="")

lines(locpoly(multipole[1:40],power[1:40],bandwidth =hn,deg=0),lwd=3,col="red")

plot(multipole[1:40],power[1:40], main="Locally polynomials with p=1",ylab="",xlab="")

lines(locpoly(multipole[1:40],power[1:40],bandwidth =hn,deg=1),lwd=3,col="red")


#################
####exmapel 6####
#################

r=function(x){sqrt(x*(1-x))*sin((2.1*pi)/(x+.05))}
curve(r,main="Doppler function")
i=seq(0,1,length=1000)

yi=r(i)+.1*rnorm(1000,0,1)
i=seq(0,1,length=1000)
plot(i,yi,main="the data")


hcv(i,yi,display="lines",ngrid=80)

















plot(i,yi,main="fitted function")
lines(ksmooth(i,yi, "normal", bandwidth=0.005),col=6,lwd=2)
lines(ksmooth(i,yi, "normal", bandwidth=0.07),col=3,lwd=2)
hcv(i,yi, display="lines",ngrid=50)  ##ostad in dastor ba moshkel hamrah hast##
a=sapply(i,hcv(yi, display="lines",ngrid=50))
plot(a)



















library(CVThresh)
?dopp