#install.packages("rngWELL")

#install.packages("randtoolbox")

library(rngWELL)

library(randtoolbox)

halt=halton(10,dim=1)

halt

write.csv(halt,"halt.csv")

> polyroot(c(40,0,0,5,20,10))
[1]  0.7747767+0.7263645i -0.7747767+1.0830293i -0.7747767-1.0830293i
[4]  0.7747767-0.7263645i -2.0000000+0.0000000i
> 

install.packages("cubature")
library(cubature)

> testFn0 <- function(x) {
+   prod(cos(x))
+ }
> 
> adaptIntegrate(testFn0, rep(0,2), rep(1,2), tol=1e-4)
$integral
[1] 0.7080734

$error
[1] 1.709434e-05

$functionEvaluations
[1] 17

$returnCode
[1] 0



> M_2_SQRTPI <- 2/sqrt(pi)
> testFn1 <- function(x) {
+   scale = 1.0
+   val = 0
+   dim = length(x)
+   val = sum (((1-x) / x)^2)
+   scale = prod(M_2_SQRTPI/x^2)
+   exp(-val) * scale
+ }
> 
> adaptIntegrate(testFn1, rep(0, 3), rep(1, 3), tol=1e-4)
$integral
[1] 1.00001

$error
[1] 9.677977e-05

$functionEvaluations
[1] 5115

$returnCode
[1] 0




> testFn2 <- function(x) {
+   ## discontinuous objective: volume of hypersphere
+   radius = as.double(0.50124145262344534123412)
+   ifelse(sum(x*x) < radius*radius, 1, 0)
+ }
> 
> adaptIntegrate(testFn2, rep(0, 2), rep(1, 2), tol=1e-4)
$integral
[1] 0.19728

$error
[1] 1.972614e-05

$functionEvaluations
[1] 166141

$returnCode
[1] 0

> adaptIntegrate(testFn3, rep(0,3), rep(1,3), tol=1e-4)
$integral
[1] 1

$error
[1] 2.220446e-16

$functionEvaluations
[1] 33

$returnCode
[1] 0


> testFn4 <- function(x) {
+   a = 0.1
+   s = sum((x-0.5)^2)
+   (M_2_SQRTPI / (2. * a))^length(x) * exp (-s / (a * a))
+ }
> 
> adaptIntegrate(testFn4, rep(0,2), rep(1,2), tol=1e-4)
$integral
[1] 1.000003

$error
[1] 9.843987e-05

$functionEvaluations
[1] 1853

$returnCode
[1] 0



> testFn5 <- function(x) {
+   a = 0.1
+   s1 = sum((x-1/3)^2)
+   s2 = sum((x-2/3)^2)
+   0.5 * (M_2_SQRTPI / (2. * a))^length(x) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
+ }
> 
> adaptIntegrate(testFn5, rep(0,3), rep(1,3), tol=1e-4)
$integral
[1] 0.9999937

$error
[1] 9.980147e-05

$functionEvaluations
[1] 59631

$returnCode
[1] 0



> testFn6 <- function(x) {
+   a = (1+sqrt(10.0))/9.0
+   prod(a/(a+1)*((a+1)/(a+x))^2)
+ }
> 
> adaptIntegrate(testFn6, rep(0,4), rep(1,4), tol=1e-4)
$integral
[1] 0.9999984

$error
[1] 9.996851e-05

$functionEvaluations
[1] 18753

$returnCode
[1] 0



> testFn7 <- function(x) {
+   n <- length(x)
+   p <- 1/n
+   (1+p)^n * prod(x^p)
+ }
> adaptIntegrate(testFn7, rep(0,3), rep(1,3), tol=1e-4)
$integral
[1] 1.000012

$error
[1] 9.966567e-05

$functionEvaluations
[1] 7887

$returnCode
[1] 0

> I.1d <- function(x) {

+   sin(4*x) *

+     x * ((x * ( x * (x*x-4) + 1) - 1))

+ }

> 

> adaptIntegrate(I.1d, -2, 2, tol=1e-7)

$integral

[1] 1.635644

$error

[1] 4.024021e-09

$functionEvaluations

[1] 105

$returnCode

[1] 0

> adaptIntegrate(I.2d, rep(-1, 2), rep(1, 2), maxEval=10000)

$integral

[1] -0.01797993

$error

[1] 7.845607e-07

$functionEvaluations

[1] 10013

$returnCode

[1] 0

> dmvnorm <- function (x, mean, sigma, log = FALSE) {
+     if (is.vector(x)) {
+         x <- matrix(x, ncol = length(x))
+     }
+     if (missing(mean)) {
+         mean <- rep(0, length = ncol(x))
+     }
+     if (missing(sigma)) {
+         sigma <- diag(ncol(x))
+     }
+     if (NCOL(x) != NCOL(sigma)) {
+         stop("x and sigma have non-conforming size")
+     }
+     if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
+         check.attributes = FALSE)) {
+         stop("sigma must be a symmetric matrix")
+     }
+     if (length(mean) != NROW(sigma)) {
+         stop("mean and sigma have non-conforming size")
+     }
+     distval <- mahalanobis(x, center = mean, cov = sigma)
+     logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
+     logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
+     if (log)
+         return(logretval)
+     exp(logretval)
+ }
> 
> m <- 3
> sigma <- diag(3)
> sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
> sigma[3,2] <- sigma[2, 3] <- 11/15
> adaptIntegrate(dmvnorm, lower=rep(-0.5, m), upper=c(1,4,2),
+                         mean=rep(0, m), sigma=sigma, log=FALSE,
+                maxEval=10000)
$integral
[1] 0.3341125

$error
[1] 4.185435e-06

$functionEvaluations
[1] 10065

$returnCode
[1] 0

g <- function(x) exp(-x)

MC.simple.est(g, 2, 4)

 [1] 0.1161926

> MC.simple.est <- function(g, a, b, n=1e4) {
+   xi <- runif(n,a,b)      # step 1
+   g.mean <- mean(g(xi))   # step 2
+   (b-a)*g.mean            # step 3
+ }
> g <- function(x) 1/log(x)
> 
> MC.simple.est(g, 2, 4)
[1] 1.922819
> integrate(g,2,4)
1.922421 with absolute error < 7.2e-14
>

peg.tab

peg.df

expand.gird()

> pag.tab <- matrix(c(762, 484, 327, 239, 468, 477), nrow=2)

> dimnames(pag.tab) <-list(Gender=c("Female","Male"),Party=c("Democrat","Independent","Republican"))

> pag.tab <- as.table(pag.tab)

> pag.tab

        Party

Gender   Democrat Independent Republican

  Female      762         327        468

  Male        484         239        477

> # Or

> pag.df <-expand.grid(Gender=c("Female","Male"),Party=c("Democrat","Independent","Republican"))

> pag.df

  Gender       Party

1 Female    Democrat

2   Male    Democrat

3 Female Independent

4   Male Independent

5 Female  Republican

6   Male  Republican

>

rm(list=ls())

mydata<-read.csv(file.choose())

mydata

fix(mydata)

install.packages("aplpack")

library(aplpack)

data =mydata

faces(data[1:30, c("Jamiyat", "M.R.S.J", "S.J.K", "M.SH","Masahat")], 

face.type = 1, scale = TRUE, 

labels = data$Species, 

plot.faces = TRUE, nrow.plot =6, 

ncol.plot = 5)

a<-2.7; b<-6.3; size<-1e4

f  <- function(x)   dbeta(x,a,b)
rg <- function(x)   runif(1,0,1)
g  <- function(x,y) 1 # i.e., dunif(x,0,1)

X <- metropolis.hastings(f,g,rg,x0=runif(1,0,1),chain.size=size)

par(mfrow=c(1,2),mar=c(2,2,1,1))
hist(X,breaks=50,col="blue",main="Metropolis-Hastings",freq=FALSE)
curve(dbeta(x,a,b),col="sienna",lwd=2,add=TRUE)
hist(rbeta(size,a,b),breaks=50,col="grey",main="Direct Sampling",freq=FALSE)
curve(dbeta(x,a,b),col="sienna",lwd=2,add=TRUE)

par("bg=blue4")
plot(table(rpois(100, 5)), type = "h", col = "red", lwd = 10,
     main = "rpois(100, lambda = 5)")

type="h"

col="red"

lwd=10

main = "rpois(100, lambda = 5)"

par(bg="gold")
attach(mtcars)
plot(wt, mpg, main="Scatterplot Example",
   xlab="Car Weight ", ylab="Miles Per Gallon ", pch=20) 

mtcars

xlab=" " ,ylab=" "

pch=" "

# Define the cars vector with 5 values

cars <- c(1, 3, 6, 4, 9)

# Graph the cars vector with all defaults

plot(cars)

# Define 2 vectors

cars <- c(1, 3, 6, 4, 9)

trucks <- c(2, 5, 4, 5, 12)

# Graph cars using a y axis that ranges from 0 to 12

plot(cars, type="o", col="blue", ylim=c(0,12))

# Graph trucks with red dashed line and square points

lines(trucks, type="o", pch=22, lty=2, col="red")

# Create a title with a red, bold/italic font

title(main="Autos", col.main="red", font.main=4)

# Define 2 vectors

cars <- c(1, 3, 6, 4, 9)

trucks <- c(2, 5, 4, 5, 12)

# Calculate range from 0 to max value of cars and trucks

g_range <- range(0, cars, trucks)

# Graph autos using y axis that ranges from 0 to max 

# value in cars or trucks vector.  Turn off axes and 

# annotations (axis labels) so we can specify them ourself

plot(cars, type="o", col="blue", ylim=g_range, 

   axes=FALSE, ann=FALSE)

# Make x axis using Mon-Fri labels

axis(1, at=1:5, lab=c("Mon","Tue","Wed","Thu","Fri"))

# Make y axis with horizontal labels that display ticks at 

# every 4 marks. 4*0:g_range[2] is equivalent to c(0,4,8,12).

axis(2, las=1, at=4*0:g_range[2])

# Create box around plot

box()

# Graph trucks with red dashed line and square points

lines(trucks, type="o", pch=22, lty=2, col="red")

# Create a title with a red, bold/italic font

title(main="Autos", col.main="red", font.main=4)

# Label the x and y axes with dark green text

title(xlab="Days", col.lab=rgb(0,0.5,0))

title(ylab="Total", col.lab=rgb(0,0.5,0))

# Create a legend at (1, g_range[2]) that is slightly smaller 

# (cex) and uses the same line colors and points used by 

# the actual plots 

legend(1, g_range[2], c("cars","trucks"), cex=0.8, 

   col=c("blue","red"), pch=21:22, lty=1:2)