1.
library(ggplot2) library(cubature)
f = function(x) {dt(x,1,0)} w = rt(10000,1,0) ggplot(w,aes(x)) + geom_density() + stat_function(fun = f, color = "red")
betaPDF <- function(x) { ifelse(0 < x & x < 1, 4*x^3, 0)} bpi = function(x){ 1/pi/(1+x^2) } oneDsampleplot <- function(oneDsample){ w<-data.frame(oneDsample) ggplot(w,aes(x)) + geom_density() + stat_function(fun = f, color = "red") }
bdtest <- runif(200000,-50,50) lb = min(bdtest[which(f(bdtest)>0)]) ub = max(bdtest[which(f(bdtest)>0)])
optimize(betaPDF,c(-1,10),maximum = TRUE)
oneDsample <- function(f, N, lb = -Inf, ub = Inf, mean = 0) { if (abs(integrate(f, lb, ub)$val - 1) > 0.001) { stop("Error: Bound is missing/wrong or the function is not a pdf. The area under the function you given should be 1") } else{ if (lb != -Inf & ub != Inf){ maxf <- optimize(f,c(lb,ub),maximum = TRUE) maxf <- maxf$objective ones = c() n = 0 while (n < N) { one <-runif(1,lb,ub) if (runif(1, 0, maxf) < f(one)){ ones = c(ones, one) n = n + 1 } } return(data.frame(x=ones)) } else { x <- runif(200000,-5000,5000) maxf <- max(f(x)) mu=x[which(f(x)==maxf)] sd = 2/maxf C = maxf/dnorm(mu,mu,sd) ones = c() n = 0 while (n < N) { one = rnorm(1, mu, sd) if (runif(1, 0, C * dnorm(one,mu,sd)) < f(one)){ ones = c(ones, one) n = n + 1 } } return(data.frame(x=ones)) } } }
f<- function(x) {dt(x,1)} f<- function(x) {dunif(x,-100,100)} f<- function(x) {ifelse(x>0, 0.0005*exp(-0.0005*x),0)} f<- function(x){ifelse(9998 < x & x < 10000, 2/39996*x, 0)} f<- function(x){ifelse(0 < x & x < 4, 2/16*x, 0)} f<- function(x) dnorm(x,-1000,10000) f<- function(x) 1/pi/(1+x^2) f<- function(x){ifelse(0 < x & x < 1, 4*x^3, 0)} f<- function(x) {ifelse(0<=x, dlnorm(x,mean=0,sdlog=1),0)} integrate(f,-Inf,Inf) a<-oneDsample(f,20000) a<-oneDsample(f,20000,-100,100) optimize(f,c(-5000,5000),maximum = TRUE) optim(jointPFF) oneDsampleplot(a) fr <- function(x) { ## Rosenbrock Banana function x1 <- x[1] x2 <- x[2] 100 * (x2 - x1 * x1)^2 + (1 - x1)^2 } (c(-1.2,1), fr)
2.
jointPFF <- function(x){ x1 = x[1] x2 = x[2] ifelse(0<x1 & x1 <1 & 0<x2 & x2<1 & 0<x1+x2 & x1+x2<1, 24*x1*x2, 0)} f <- function(x){ x1 = x[1] x2 = x[2] dnorm(x1,0,1) * dnorm(x2, 2, 3) } j2 <- function(x){ x1 = x[1] x2 = x[2] 1/pi/(1+x1^2) * 0.05*exp(-0.05*x2) } ggplot(twoDsample(f, N=5000), aes(x, y)) + geom_density_2d() mmmv = mvrnorm(100,c(0,0),matrix(c(1,0,0,1),2,2)) mv = data.frame(x = mmmv[,1], y = mmmv[,2]) ggplot(mv, aes(x,y)) + geom_density_2d() optim(c(0.3,0.2),j1, control = list(fnscale = -1)) nlm(jointPFF,c(0,9)) adaptIntegrate(j2, c(-1000,0), c(1000,1000))
twoDsample <- function(f, N, lbx, ubx, lby, uby) { library(cubature) if (missing(lbx)){ lbx = -5000 } if (missing(lby)){ lby = -5000 } if (missing(ubx)){ ubx = 5000 } if (missing(uby)){ uby = 5000 } #if (abs(adaptIntegrate(f, c(lbx, lby), c(ubx, uby), maxEval=10000)$integral - 1) > 0.001) { #stop("Error: Bound is missing/wrong or the function is not a pdf. The area under the function you given should be 1") #} #else{ # op = optim(c((ubx+lbx)/2,(uby+lby)/2),f, control = list(fnscale = -1)) # twos = c() # n = 0 # while (n < N) { # two <- c(runif(1,lbx,ubx),runif(1,lby,uby)) # if (runif(1, 0, maxf) < f(two)){ # twos = c(twos, two) # n = n + 1 #} #} #return(data.frame(x=twos[c(seq(1,length(twos)-1,2))],y=twos[c(seq(2,length(twos),2))])) #} #else{ op = optim(c((ubx+lbx)/2,(uby+lby)/2), j2, control = list(fnscale = -1)) maxf = op$value mu = c(op$par) sd = 2/maxf C = maxf/dmvnorm(c(mu[1],mu[2]),c(mu[1],mu[2]),c(sd,sd)) twos = c() n = 0 while (n < N) { two = mvrnorm(1, mu, matrix(c(sd,0,0,sd),2,2)) if (runif(1, 0, C * dmvnorm(two,mu,c(sd,sd))) < f(two)){ twos = c(twos, two) n = n + 1 } } return(data.frame(x=twos[c(seq(1,length(twos)-1,2))],y=twos[c(seq(2,length(twos),2))])) #} }
box = c() maxf = 0 for (i in 1:1000000){ xy = c(runif(1,-1000,1000),runif(1,0,1000)) if (maxf < j2(xy)){ maxf = j2(xy) box = xy } } maxf xy Sigma
maxf <- max(replicate(200000,j2(c(runif(1,-50,50),runif(1,0,5000))))) xy <- c(runif(200000,-5000,5000), y = runif(200000,0,5000)) if (jointPFF(xy) == maxfxy){ x = x } x
x1 = oneDsample(f = betaPDF, N=100, lb = 0, ub = 1) x2 = twoDsample(fj = jointPFF, N=100, lbx=0, ubx=1, lby=0, uby=1) g1 = function(x){x^2} g2 = function(x,y){x^2+y} twoEvalue(g1,x1) twoEvalue(g2,x2)
mean(x1<0.5) gb = function(x){ x<0.5 } mean(gb(x1)) c1 = function(x){x<0.4} c2 = function(x,y){ x<0.5 & y<0.9 } Gprob(c1,x1) Gprob(c2,x2)
mmmv = mvrnorm(10000,c(55,55),matrix(c(7,0,0,7),2,2)) mv = data.frame(x = mmmv[,1], y = mmmv[,2]) #plot(mv) ggplot(mv, aes(x,y)) + geom_density_2d() + geom_abline_2d() ggplot(mv, aes(x,y)) + geom_density2d()
mmmv = mvrnorm(10000,c(55,55),matrix(c(1,0,0,1),2,2)) mv = data.frame(x = mmmv[,1], y = mmmv[,2]) dmv = dmvnrom(mv,c(55,55),c(1,1)) max(dmv) ggplot(mv, aes(x)) + geom_density()
dmvnorm = function(x,mu,sig){ x1 = x[1] x2 = x[2] mu1 = mu[1] mu2 = mu[2] sig1 = sig[1] sig2 = sig[2] exp(-1/2*((x1-mu1)^2/sig1^2 - 2*(x1-mu1)*(x2-mu2)/sig1/sig2 + (x2-mu2)^2/sig2^2))/(2*pi*sig1*sig2) }
f <- function(x){ x1 = x[1] x2 = x[2] ifelse(0<x1 & x1 <1 & 0<x2 & x2<1 & 0<x1+x2 & x1+x2<1, 24*x1*x2, 0)} if (f(c(-50,-50)) == 0 & f(c(50,50)) == 0){ flag = 0 lbx = -50 lby = -50 ubx = 50 uby = 50 while (flag == 0){ min = optim(c((ubx+lbx)/2,(uby+lby)/2), f) if (min$value == 0){ lbx = lbx/2 lby = lby/2 ubx = ubx/2 uby = uby/2 } if (min$value > 0){ lbx = lbx*1.1 lby = lby*1.1 ubx = ubx*1.1 uby = uby*1.1 } if (f(c(lbx-0.01,lby-0.01)) == 0 & f(c(lbx,lby)) > 0 & f(c(ubx+0.01,uby+0.01)) == 0 & f(c(ubx,uby)) > 0){ flag = 1 } } }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.