# In Kitetsu3rd/Project1:

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)
}

f <- function(x){
x1 = x[1]
x2 = x[2]
1/pi/(1+x1^2) * 0.05*exp(-0.05*x2)
}

f <- function(x){
x1 = x[1]
x2 = x[2]
1/pi/(1+x1^2) * dt(x2,1)
}
ggplot(twoDsample(f, N=5000), aes(x, y)) +  geom_density_2d()
ggplot(twoDsample(f, N=5000,lby = 0,lbx=0,ubx=1,uby=1), aes(x, y)) +  geom_density_2d()
```
```ggplot(twoDsample(j2, N=5000,lby = 0), 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.1,0.1),jointPFF, control = list(fnscale = -1))
nlm(jointPFF,c(0,9))
adaptIntegrate(f, c(-5000,-5000), c(5000,5000), maxEval = 100000)
```
```twoDsample <- function(f, N, lbx=-5000, ubx=5000, lby=-5000, uby=5000) {
library(cubature)
library(MASS)
if (abs(adaptIntegrate(f, c(lbx, lby), c(ubx, uby), maxEval=100000)\$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))
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)
}
```
```optim(c(-2,3), j2, control = list(fnscale = -1))
```

Kitetsu3rd/Project1 documentation built on May 30, 2019, 3:49 p.m.