################################################################################
#	        R code to construct histogram of PM2.5 data           #
################################################################################
Maehongson <- c(25,21,19,18,18,20,16,23,25,26,20,18,17,21,27,28,24,28,27,28,25,21,26,30,36,38,40,42,41,35,40,41,47,42,41,34,30,31,37,39,54,64,78,69,57,70,68,13,9,15,21,27,28,34,36,48,59,66,75,78,79,92,97,103,104,91,91,111,118,142,216,134,53,59,62,60,98,116,104,148,137,144,255,309,244,237,209,263,321,256,233,144,109,147,176,204,203,147,117,108,124,119,104,106,166,214,138,61,59,56,53,50,38,34,34,49,39,36,50,49,34,23,23,22,28,29,31,26,13,18,8,8,5,11,12,13,17,18,19,19,21,19,19,15,17,17,16,17,13,10,11,8.3,11.6,13.3,10.3,6.2,5.5,5.8,7.1,14.1,13.5,15.2,10,10.1,10.2,10.2,9,5.5,3.7,3.7,3.9,3.6,3.3,4.2,4.2,4.4,4.4,5,5.5,4.8,7.1)
h <- hist(Maehongson,xlab=" PM2.5 levels",xlim=c(0,350),ylim=c(0,60), main="(A) Mae Hong Son province")
xfit <- seq(min(Maehongson),max(Maehongson),length=50)
yfit <- dnorm(xfit,mean=mean(Maehongson),sd=sd(Maehongson))
yfit <- yfit*diff(h$mids[1:2])*length(Maehongson)
lines(xfit,yfit,col="blue",lwd=2)
Lampang <- c(29,30,28,30,27,21,19,33,39,32,37,33,30,27,32,41,43,39,41,40,44,50,49,57,60,59,71,62,59,67,88,87,80,82,68,52,40,39,65,106,108,113,125,158,211,153,42,23,19,27,57,67,74,72,88,92,59,96,143,155,150,138,132,119,100,84,90,97,101,99,103,64,49,48,53,41,51,34,32,39,42,51,62,129,119,99,121,94,103,148,148,86,64,92,172,185,182,113,80,72,96,133,130,100,122,119,82,89,91,93,91,60,55,37,61,48,30,38,44,34,36,63,32,37,48,52,51,46,17,19,21,20,16,17,22,23,23,35,33,31,36,27,24,23,29,31,28,34,25,22,18,18.6,17.4,20.3,20.6,16.1,16.8,16.7,14.9,16.5,19.2,20.4,21.6,18.7,16.91,8.2,16.3,16,14.8,15.6,14.8,13.5,13.3,14.7,12.8,13,12.5,12.7,14.4,13,12.7)
h <- hist(Lampang,xlab=" PM2.5 levels",xlim=c(0,350),ylim=c(0,60), main="(B) Lampang province")
xfit <- seq(min(Lampang),max(Lampang),length=50)
yfit <- dnorm(xfit,mean=mean(Lampang),sd=sd(Lampang))
yfit <- yfit*diff(h$mids[1:2])*length(Lampang)
lines(xfit,yfit,col="blue",lwd=2)