stuff

library(ggplot2)

# Set parameters
odin.model <- function(x)
  sqrt(x+0.3)
mean.1 <- 0.75
mean.2 <- 1.25
pm.sd <- 0.4
odin.sd <- 0.1
ss.1 <- 100 # sample size
ss.2 <- 100

# Generate the PM concnetrations
pm.1 <- rnorm(ss.1,mean.1,pm.sd)
pm.2 <- rnorm(ss.2,mean.2,pm.sd)

# Generate the odin measurements
odin.1 <- odin.model(pm.1) + rnorm(ss.1,0,odin.sd)
## Warning in sqrt(x + 0.3): NaNs produced
odin.2 <- odin.model(pm.2) + rnorm(ss.2,0,odin.sd)

# Create dataframe
data.1 <- data.frame(pm=pm.1,odin=odin.1,data=rep('sample.1',ss.1))
data.2 <- data.frame(pm=pm.2,odin=odin.2,data=rep('sample.2',ss.2))
data <- rbind(data.1,data.2)

# Generate linear models
lin.model.1 <- lm(odin~pm,data=data.1)
lin.model.2 <- lm(odin~pm,data=data.2)
line.1 <- data.frame(pm=c(mean.1-2*pm.sd,mean.1+2*pm.sd),odin=NA,data=rep('sample.1',2))
line.1[,'odin'] <- predict(lin.model.1,newdata=line.1)
line.2 <- data.frame(pm=c(mean.2-2*pm.sd,mean.2+2*pm.sd),odin=NA,data=rep('sample.2',2))
line.2[,'odin'] <- predict(lin.model.2,newdata=line.2)
lines <- rbind(line.1,line.2)

# Distribution visualizations
pm.range <- seq(-0.5,2.5,length.out=100) 
dist.1 <- data.frame(pm=pm.range,odin=dnorm(pm.range,mean.1,pm.sd),
                     data=rep('sample.1',100))
dist.2 <- data.frame(pm=pm.range,odin=dnorm(pm.range,mean.2,pm.sd),
                     data=rep('sample.2',100))
dists <- rbind(dist.1,dist.2)

# ODIN model visualization
odin.curve <- data.frame(pm=pm.range,odin=odin.model(pm.range))
## Warning in sqrt(x + 0.3): NaNs produced
# Plot
ggplot() +
  geom_point(data=data,aes(pm,odin,colour=data),alpha=0.5) +
  geom_line(data=dists,aes(pm,odin,colour=data)) +
  geom_line(data=lines,aes(pm,odin,colour=data),size=1) +
  geom_line(data=odin.curve,aes(pm,odin),linetype="dashed",size=0.75,colour="black")
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 7 rows containing missing values (geom_path).

Weights according to density function

# Calculate weights
weights.1 <- dnorm(pm.1,mean.2,pm.sd)
weights.2 <- dnorm(pm.2,mean.1,pm.sd)

# Add weights to data
data[,'weight'] <- c(weights.1,weights.2)

# Geneate linear models
lin.model.1 <- lm(odin~pm,data=data.1,weights=weights.1)
lin.model.2 <- lm(odin~pm,data=data.2,weights=weights.2)
line.1 <- data.frame(pm=c(mean.1-2*pm.sd,mean.1+2*pm.sd),odin=NA,data=rep('sample.1',2))
line.1[,'odin'] <- predict(lin.model.1,newdata=line.1)
line.2 <- data.frame(pm=c(mean.2-2*pm.sd,mean.2+2*pm.sd),odin=NA,data=rep('sample.2',2))
line.2[,'odin'] <- predict(lin.model.2,newdata=line.2)
lines <- rbind(line.1,line.2)

# Distribution visualizations
pm.range <- seq(-0.5,2.5,length.out=100) 
dist.1 <- data.frame(pm=pm.range,odin=
                    dnorm(pm.range,mean.1,pm.sd)*dnorm(pm.range,mean.2,pm.sd),
                     data=rep('sample.1',100))
dist.2 <- data.frame(pm=pm.range,odin=
                    dnorm(pm.range,mean.1,pm.sd)*dnorm(pm.range,mean.2,pm.sd),
                     data=rep('sample.2',100))
dists <- rbind(dist.1,dist.2)

# ODIN model visualization
odin.curve <- data.frame(pm=pm.range,odin=odin.model(pm.range))
## Warning in sqrt(x + 0.3): NaNs produced
# Plot
ggplot() +
  geom_point(data=data,aes(pm,odin,colour=data,size=weight),alpha=0.25) +
  scale_size_continuous(range=c(0.1,2.5)) +
  geom_line(data=dists,aes(pm,odin,colour=data)) +
  geom_line(data=lines,aes(pm,odin,colour=data),size=1) +
  geom_line(data=odin.curve,aes(pm,odin),linetype="dashed",size=0.75,colour="black")
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 7 rows containing missing values (geom_path).

Weights according to binning

bin.breaks <- seq(-2,2.5,length.out=20)
pm.bin.1 <- cut(pm.1,breaks=bin.breaks)
pm.bin.2 <- cut(pm.2,breaks=bin.breaks)

# The bumber of occurances in each bin
pm.num.bin.1 <- summary(pm.bin.1)
pm.num.bin.2 <- summary(pm.bin.2)

# The density of pm concentratinos from the two samples (binned)
pm.dist.1 <- summary(pm.bin.1)/sum(pm.num.bin.1)
pm.dist.2 <- summary(pm.bin.2)/sum(pm.num.bin.2)

# The desired distribution (binned)
#pm.dist.out <- data.frame(pm=pm.dist.1)# * pm.dist.2)

# Calculate weights
weights.2 <- pm.dist.1[pm.bin.2]
weights.1 <- pm.dist.2[pm.bin.1]

# Add weights to data
data[,'weight'] <- c(weights.1,weights.2)

# Geneate linear models
lin.model.1 <- lm(odin~pm,data=data.1,weights=weights.1)
lin.model.2 <- lm(odin~pm,data=data.2,weights=weights.2)
line.1 <- data.frame(pm=c(mean.1-2*pm.sd,mean.1+2*pm.sd),odin=NA,data=rep('sample.1',2))
line.1[,'odin'] <- predict(lin.model.1,newdata=line.1)
line.2 <- data.frame(pm=c(mean.2-2*pm.sd,mean.2+2*pm.sd),odin=NA,data=rep('sample.2',2))
line.2[,'odin'] <- predict(lin.model.2,newdata=line.2)
lines <- rbind(line.1,line.2)

# Distribution visualizations
pm.range <- seq(-0.5,2.5,length.out=100) 
dist.1 <- data.frame(pm=pm.range,odin=
                    dnorm(pm.range,mean.1,pm.sd)*dnorm(pm.range,mean.2,pm.sd),
                     data=rep('sample.1',100))
dist.2 <- data.frame(pm=pm.range,odin=
                    dnorm(pm.range,mean.1,pm.sd)*dnorm(pm.range,mean.2,pm.sd),
                     data=rep('sample.2',100))
dists <- rbind(dist.1,dist.2)

# ODIN model visualization
odin.curve <- data.frame(pm=pm.range,odin=odin.model(pm.range))
## Warning in sqrt(x + 0.3): NaNs produced
# Plot
ggplot() +
  geom_point(data=data,aes(pm,odin,colour=data,size=weight),alpha=0.25) +
  scale_size_continuous(range=c(0.1,2.5)) +
  geom_line(data=dists,aes(pm,odin,colour=data)) +
  geom_line(data=lines,aes(pm,odin,colour=data),size=1) +
  geom_line(data=odin.curve,aes(pm,odin),linetype="dashed",size=0.75,colour="black") #+
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 7 rows containing missing values (geom_path).

  #geom_histogram(breaks=bin.breaks,data=pm.dist.out,aes(pm))