# rough plot but with full data:
n=10^5
xin=runif(n)
yin=runif(n)
f=runif(n)*5
x=xin[f<2*xin+3*yin^2]
y=yin[f<2*xin+3*yin^2]
xvec=seq(0,1,l=400)
ertol=0.02
muy=numeric()
for (k in 1:length(xvec)){
    muy[k]=mean(y[x>xvec[k]-ertol & x<xvec[k]+ertol])
}
plot(x,y,type="p",pch=20,col="grey",cex=0.8)
lines(xvec,muy,type="p",pch=20,col="black",cex=1.2)
lines(xvec,(xvec+3/4)/(2*xvec+1),type="l",col="red",lwd=2)


# well-fitted plot:
# may take a long time to run w/ n=10^6
n=10^6
xin=runif(n)
yin=runif(n)
f=runif(n)*5
x=xin[f<2*xin+3*yin^2]
y=yin[f<2*xin+3*yin^2]
xvec=seq(0,1,l=400)
ertol=0.02
muy=numeric()
for (k in 1:length(xvec)){
    muy[k]=mean(y[x>xvec[k]-ertol & x<xvec[k]+ertol])
}
plot(xvec,muy,type="p",pch=20,col="black",cex=1.2,xlab="x",ylab="mu_y|x = E[Y|X=x]")
lines(xvec,(xvec+3/4)/(2*xvec+1),type="l",col="red",lwd=2)